%% Assess Model Stability Using Rolling Window Analysis
% This example shows how to use a rolling window to check whether the parameters of a time-series model are time invariant. This example analyzes two time series:
% Time-series 1: simulated data from a known, time-invariant model
% Time-series 2: simulated data from a known, time-varying model
% Completely specify this AR(1) model for Time-series 1:
% y_t = 0.6 y_{t-1} + \varepsilon_t
% where \varepsilon_t is Gaussian with mean 0 and variance 1. Completely specify this time-varying model for Time-series 2:
% \begin{array}{*{20}{c}}
% {{y_t} = 0.2{y_{t - 1}} + {\varepsilon _t}{\rm{; }}t = 1,...,100}\\
% {{y_t} = 0.75{y_{t - 1}} + {\varepsilon _t}{\rm{; }}t = 101,...,150}\\
% {{y_t} =  - 0.5{y_{t - 1}} + {\varepsilon _t}{\rm{; }}t = 151,...,200},
% \end{array}
% where \varepsilon_t is Gaussian with mean 0 and variance 1.

Mdl1 = arima('AR',0.6,'Constant',0,'Variance',1);

Mdl2 = cell(3,1); % Preallocate
ARMdl2 = [0.2 0.75 -0.5];
for j = 1:3
    Mdl2{j} = arima('AR',ARMdl2(j),'Constant',0,'Variance',1);
end

% Mdl1 is an arima model objects. You can access its properties using dot notation.  Mdl2 is a cell array of arima model objects. You can you cell indexing and dot notation to access properties of the models within Mdl2. For example, to access the AR parameter value of the third model in Mdl3, enter Mdl2{3}.AR.
%
% Simulate T = 200 periods of data from Mdl1 and Mdl2. Use a presample response of 0 for both series.

rng(1); % For reproducibility
T = 200;
y1 = simulate(Mdl1,T,'Y0',0);
timeMdl2 = [100 50 50]; % Number of observations per model in Mdl2
y2 = 0;
for k = 1:numel(Mdl2)
    y2 = [y2; simulate(Mdl2{k},timeMdl2(k),'Y0',y2(end))];
end

Y = [y1 y2(2:end)];

%
% Specify empty AR(1) models for the estimation of Mdl1, Mdl2, and Mdl3. Estimate all three models using the respective data sets and a rolling window size of 40 periods. Also, use a rolling window increment of one period. Store the autoregressive parameters and estimated innovations variance.

Mdl = arima('ARLags',1,'Constant',0);

m = 100;       % Rolling window size
N = T - m + 1; % Number of rolling windows

EstParams = cell(2,1);     % Preallocate for estimates
EstParamsMat = zeros(N,2);
EstParamsSE = cell(2,1);
EstParamsSEMat = zeros(N,2);

for j = 1:2
    for k = 1:N
        idxRW =  k:(m + k - 1); % In-sample indices
        [EstMdl,EstParamCov] = estimate(Mdl,Y(idxRW,j),'Display','off');
        EstParamsMat(k,:) = [EstMdl.AR{1} EstMdl.Variance];
        EstParamsSEMat(k,:) = sqrt([EstParamCov(2,2) EstParamCov(3,3)]);
    end
    EstParams{j} = EstParamsMat;
    EstParamsSE{j} = EstParamsSEMat;
end

%
% Plot the estimates and their point-wise confidence intervals over the rolling window index.

titleMdls = {'Mdl1','Mdl2'};

for j = 1:2
    figure;
    subplot(2,1,1);
    Estimates = EstParams{j};
    SEs = EstParamsSE{j};
    plot(Estimates(:,1),'LineWidth',2);
    hold on;
    plot(Estimates(:,1) + 2*SEs(:,1),'r:','LineWidth',2);
    plot(Estimates(:,1) - 2*SEs(:,1),'r:','LineWidth',2);
    title(sprintf('%s - AR at Lag 1 Estimate',titleMdls{j}));
    xlabel 'Rolling window index';
    axis tight;
    hold off;
    subplot(2,1,2);
    plot(Estimates(:,2),'LineWidth',2);
    hold on;
    plot(Estimates (:,2) + 2*SEs(:,2),'r:','LineWidth',2);
    plot(Estimates(:,2) - 2*SEs(:,2),'r:','LineWidth',2);
    title(sprintf('%s - Variance Estimate',titleMdls{j}));
    xlabel 'Rolling window index';
    axis tight;
    hold off;
end

% For Mdl1, the AR estimate does not vary much from 0.6, and the estimates are not significantly different from one another (pair-wise). Similar results occur for the variance of Mdl1. The AR estimate of Mdl2 grows, and then falls, which indicates time dependence. Also, based on the confidence intervals, there is evidence that some estimates differ from others. Though the variance did not change during simulation, there seems to be heteroscedasticity possibly induced by the instability of the model.
% Copyright 2015 The MathWorks, Inc.