%% Simulate an AR Error Model
% This example shows how to simulate sample paths from a regression model with AR errors without specifying presample disturbances.
% Specify the regression model with AR(2) errors:
% \begin{array}{l}{y_t} = 2 + {X_t}\left[ \begin{array}{l} - 2\\1.5\end{array} \right] + {u_t}\\{u_t} = 0.75{u_{t - 1}} - 0.5{u_{t - 2}} + {\varepsilon _t},\end{array}
% where \varepsilon_t is Gaussian with mean 0 and variance 1.

Beta = [-2; 1.5];
Intercept = 2;
a1 = 0.75;
a2 = -0.5;
Variance = 1;
Mdl = regARIMA('AR',{a1, a2},'Intercept',Intercept,...
    'Beta',Beta,'Variance',Variance);

%
% Generate two length T = 50 predictor series by random selection from the standard Gaussian distribution.

T = 50;
rng(1); % For reproducibility
X = randn(T,2);

% The software treats the predictors as nonstochastic series.
%
% Generate and plot one sample path of responses from Mdl.

rng(2);
ySim = simulate(Mdl,T,'X',X);

figure
plot(ySim)
title('{\bf Simulated Response Series}')

% simulate requires P = 2 presample unconditional disturbances (u_t) to initialize the error series. Without them, as in this case, simulate sets the necessary presample unconditional disturbances to 0.
%
% Alternatively, filter a random innovation series through Mdl using filter.

rng(2);
e = randn(T,1);
yFilter = filter(Mdl,e,'X',X);

figure
plot(yFilter)
title('{\bf Simulated Response Series Using Filtered Innovations}')

% The plots suggest that the simulated responses and the responses generated from the filtered innovations are equivalent.
%
% Simulate 1000 response paths from Mdl. Assess transient effects by plotting the unconditional disturbance (U) variances across the simulated paths at each period.

numPaths = 1000;
[Y,~,U] = simulate(Mdl,T,'NumPaths',numPaths,'X',X);

figure
h1 = plot(Y,'Color',[.85,.85,.85]);
title('{\bf 1000 Simulated Response Paths}')
hold on
h2 = plot(1:T,Intercept+X*Beta,'k--','LineWidth',2);
legend([h1(1),h2],'Simulated Path','Mean')
hold off

figure
h1 = plot(var(U,0,2),'r','LineWidth',2);
hold on
theoVarFix = ((1-a2)*Variance)/((1+a2)*((1-a2)^2-a1^2));
h2 = plot([1 T],[theoVarFix theoVarFix],'k--','LineWidth',2);
title('{\bf Unconditional Disturbance Variance}')
legend([h1,h2],'Simulation Variance','Theoretical Variance')
hold off

% The simulated response paths follow their theoretical mean, c + X\beta, which is not constant over time (and might look nonstationary).
% The variance of the process is not constant, but levels off at the theoretical variance by the 10th period. The theoretical variance of the AR(2) error model is
% \frac{{\left( {1 - {a_2}} \right)\sigma _\varepsilon ^2}}{{\left( {1 + {a_2}} \right)\left[ {{{\left( {1 - {a_2}} \right)}^2} - a_1^2} \right]}} = \frac{{(1 + 0.5)}}{{(1 - 0.5)\left[ {{{(1 + 0.5)}^2} - {{0.75}^2}} \right]}} = 1.78
%
% You can reduce transient effects is by partitioning the simulated data into a burn-in portion and a portion for analysis. Do not use the burn-in portion for analysis. Include enough periods in the burn-in portion to overcome the transient effects.

burnIn = 1:10;
notBurnIn = burnIn(end)+1:T;
Y = Y(notBurnIn,:);
X = X(notBurnIn,:);
U = U(notBurnIn,:);
figure
h1 = plot(notBurnIn,Y,'Color',[.85,.85,.85]);
hold on
h2 = plot(notBurnIn,Intercept+X*Beta,'k--','LineWidth',2);
title('{\bf 1000 Simulated Response Paths for Analysis}')
legend([h1(1),h2],'Simulated Path','Mean')
hold off

figure
h1 = plot(notBurnIn,var(U,0,2),'r','LineWidth',2);
hold on
h2 = plot([notBurnIn(1) notBurnIn(end)],...
    [theoVarFix theoVarFix],'k--','LineWidth',2);
title('{\bf Converged Unconditional Disturbance Variance}')
legend([h1,h2],'Simulation Variance','Theoretical Variance')
hold off

% Unconditional disturbance simulation variances fluctuate around the theoretical variance due to Monte Carlo sampling error. Be aware that the exclusion of the burn-in sample from analysis reduces the effective sample size.
% Copyright 2015 The MathWorks, Inc.