%% Simulate a Regression Model with Nonstationary Errors
% This example shows how to simulate responses from a regression model with ARIMA unconditional disturbances, assuming that the predictors are white noise sequences.
% Specify the regression model with ARIMA errors:
% y_t = 3 + X_t\left[\matrix{2 \cr -1.5}\right] + u_t
% \Delta u_t = 0.5\Delta u_{t-1} + \varepsilon_t + 1.4\varepsilon_{t-1} + 0.8\varepsilon_{t-2},
% where the innovations are Gaussian with variance 1.

T = 150; % Sample size
Mdl = regARIMA('MA',{1.4,0.8},'AR',0.5,'Intercept',3,...
    'Variance',1,'Beta',[2;-1.5],'D',1);

%
% Simulate two Gaussian predictor series with mean 0 and variance 1.

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

%
% Simulate and plot the response series.

y = simulate(Mdl,T,'X',X);

figure;
plot(y);
title 'Simulated Responses';
axis tight;

%
% Regress y onto X. Plot the residuals, and test them for a unit root.

RegMdl = fitlm(X,y);

figure;
subplot(2,1,1);
plotResiduals(RegMdl,'caseorder');
subplot(2,1,2);
plotResiduals(RegMdl,'lagged');

h = adftest(RegMdl.Residuals.Raw)

% The residual plots indicate that they are autocorrelated and possibly nonstationary (as constructed).  h = 0 indicates that there is insufficient evidence to suggest that the residual series is not a unit root process.
%
% Treat the nonstationary unconditional disturbances by transforming the data appropriately. In this case, difference the responses and predictors. Reestimate the regression model using the transformed responses, and plot the residuals.

dY = diff(y);
dX = diff(X);
dRegMdl = fitlm(dX,dY);

figure;
subplot(2,1,1);
plotResiduals(dRegMdl,'caseorder','LineStyle','-');
subplot(2,1,2);
plotResiduals(dRegMdl,'lagged');

h = adftest(dRegMdl.Residuals.Raw)

% The residual plots indicate that they are still autocorrelated, but stationary. h = 1 indicates that there is enough evidence to suggest that the residual series is not a unit root process.
% Once the residuals appear stationary, you can determine the appropriate number of lags for the error model using Box and Jenkins methodology. Then, use regARIMA to completely model the regression model with ARIMA errors.
% Copyright 2015 The MathWorks, Inc.