Simulate Regression Models with ARMA Errors

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} [\begin{array}{l} - 2 \\ 1.5 \end{array}] + u_{t} \\ u_{t} = 0.75 u_{t - 1} - 0.5 u_{t - 2} + ε_{t}, \end{array}$

where $ε_{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}')

Figure contains an axes object. The axes object with title blank Simulated blank Response blank Series contains an object of type line.

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}')

Figure contains an axes object. The axes object with title blank Simulated blank Response blank Series blank Using blank Filtered blank Innovations contains an object of type line.

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 contains an axes object. The axes object with title blank 1000 blank Simulated blank Response blank Paths contains 1001 objects of type line. These objects represent Simulated Path, Mean.

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

Figure contains an axes object. The axes object with title blank Unconditional blank Disturbance blank Variance contains 2 objects of type line. These objects represent Simulation Variance, Theoretical Variance.

The simulated response paths follow their theoretical mean, $c + X β$ , 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{(1 - a_{2}) σ_{ε}^{2}}{(1 + a_{2}) [{(1 - a_{2})}^{2} - a_{1}^{2}]} = \frac{(1 + 0.5)}{(1 - 0.5) [{(1 + 0.5)}^{2} - {0.75}^{2}]} = 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 contains an axes object. The axes object with title blank 1000 blank Simulated blank Response blank Paths blank for blank Analysis contains 1001 objects of type line. These objects represent Simulated Path, Mean.

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

Figure contains an axes object. The axes object with title blank Converged blank Unconditional blank Disturbance blank Variance contains 2 objects of type line. These objects represent Simulation Variance, Theoretical Variance.

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.

Simulate an MA Error Model

Open Live Script

This example shows how to simulate responses from a regression model with MA errors without specifying a presample.

Specify the regression model with MA(8) errors:

$\begin{array}{l} y_{t} = 2 + X_{t} [\begin{array}{l} - 2 \\ 1.5 \end{array}] + u_{t} \\ u_{t} = ε_{t} + 0.4 ε_{t - 1} - 0.3 ε_{t - 4} + 0.2 ε_{t - 8}, \end{array}$

where $ε_{t}$ is Gaussian with mean 0 and variance 0.5.

Beta = [-2; 1.5];
Intercept = 2;
b1 = 0.4;
b4 = -0.3;
b8 = 0.2;
Variance = 0.5;
Mdl = regARIMA('MA',{b1, b4, b8},'MALags',[1 4 8],...
    'Intercept',Intercept,'Beta',Beta,'Variance',Variance);

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

T = 100;
rng(4); % 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(5);
ySim = simulate(Mdl,T,'X',X);

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

Figure contains an axes object. The axes object with title blank Simulated blank Response blank Series contains an object of type line.

simulate requires Q = 8 presample innovations ( $ε_{t}$ ) to initialize the error series. Without them, as in this case, simulate sets the necessary presample innovations to 0.

Alternatively, use filter to filter a random innovation series through Mdl.

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

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

Figure contains an axes object. The axes object with title blank Simulated blank Response blank Series blank Using blank Filtered blank Innovations contains an object of type line.

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 contains an axes object. The axes object with title blank 1000 blank Simulated blank Response blank Paths contains 1001 objects of type line. These objects represent Simulated Path, Mean.

figure
h1 = plot(var(U,0,2),'r','LineWidth',2);
hold on
theoVarFix = (1+b1^2+b4^2+b8^2)*Variance;
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 paths follow their theoretical mean, $c + X β$ , 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 15th period. The theoretical variance of the MA(8) error model is

$(1 + b_{1}^{2} + b_{4}^{2} + b_{8}^{2}) σ_{ε}^{2} = (1 + {0.4}^{2} + {(- 0.3)}^{2} + {0.2}^{2}) 0.5 = 0.645 .$

burnIn = 1:15;
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')
axis tight
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')
axis tight
hold off

Simulate an ARMA Error Model

Open Live Script

This example shows how to simulate responses from a regression model with ARMA errors without specifying a presample.

Specify the regression model with ARMA(2,1) errors:

$\begin{array}{l} y_{t} = 2 + X_{t} [\begin{array}{l} - 2 \\ 1.5 \end{array}] + u_{t} \\ u_{t} = 0.9 u_{t - 1} - 0.1 u_{t - 2} + ε_{t} + 0.5 ε_{t - 1}, \end{array}$

where $ε_{t}$ is distributed with 15 degrees of freedom and variance 1.

Beta = [-2; 1.5];
Intercept = 2;
a1 = 0.9;
a2 = -0.1;
b1 = 0.5;
Variance = 1;
Distribution = struct('Name','t','DoF',15);
Mdl = regARIMA('AR',{a1, a2},'MA',b1,...
    'Distribution',Distribution,'Intercept',Intercept,...
    'Beta',Beta,'Variance',Variance);

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

T = 50;
rng(6); % 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(7);
ySim = simulate(Mdl,T,'X',X);

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

Figure contains an axes object. The axes object with title blank Simulated blank Response blank Series contains an object of type line.

simulate requires:

P = 2 presample unconditional disturbances to initialize the autoregressive component of the error series.
Q = 1 presample innovations to initialize the moving average component of the error series.

Without them, as in this case, simulate sets the necessary presample errors to 0.

Alternatively, use filter to filter a random innovation series through Mdl.

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

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

Figure contains an axes object. The axes object with title blank Simulated blank Response blank Series blank Using blank Filtered blank Innovations contains an object of type line.

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 contains an axes object. The axes object with title blank 1000 blank Simulated blank Response blank Paths contains 1001 objects of type line. These objects represent Simulated Path, Mean.

figure
h1 = plot(var(U,0,2),'r','LineWidth',2);
hold on
theoVarFix = Variance*(a1*b1*(1+a2)+(1-a2)*(1+a1*b1+b1^2))/...
    ((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',...
    'Location','Best')
hold off

The simulated paths follow their theoretical mean, $c + X β$ , 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 ARMA(2,1) error model is:

$\frac{σ_{ε}^{2} [a_{1} b_{1} (1 + a_{2}) + (1 - a_{2}) (1 + a_{1} b_{1} + b_{1}^{2})]}{{(1 + a_{2})}^{2} [{(1 - a_{2})}^{2} - a_{1}^{2}]} = \frac{[0.9 (0.5) (1 - 0.1) + (1 + 0.1) (1 + 0.9 (0.5) + {0.5}^{2})]}{{(1 - 0.1)}^{2} [{(1 + 0.1)}^{2} - {0.9}^{2}]} = 6 . 32 .$

You can reduce transient effects 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')
axis tight
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')
axis tight
hold off