%% Simulate a Regression Model with Nonstationary Exponential Errors
% This example shows how to simulate responses from a regression model with nonstationary, exponential, unconditional disturbances. Assume that the predictors are white noise sequences.
% Specify the following ARIMA error model:
% \Delta u_t = 0.9 \Delta u_{t-1} + \varepsilon_t,
% where the innovations are Gaussian with mean 0 and variance 0.05.

T = 50; % Sample size
MdlU = arima('AR',0.9,'Variance',0.05,'D',1,'Constant',0);

%
% Simulate unconditional disturbances. Exponentiate the simulated errors.

rng(10); % For reproducibility
u = simulate(MdlU,T,'Y0',[0.5:1.5]');
expU = exp(u);

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

X = randn(T,2);

%
% Generate responses from the regression model with time series errors:
% y_t = 3 + X_t\left[\matrix{2\cr-1.5}\right]+e^{u_t}.

Beta = [2;-1.5];
Intercept = 3;
y = Intercept + X*Beta + expU;

%
% Plot the responses.

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

% The response series seems to grow exponentially (as constructed).
%
% Regress y onto X. Plot the residuals.

RegMdl1 = fitlm(X,y);

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

% The residuals seem to grow exponentially, and seem autocorrelated (as constructed).
%
% Treat the nonstationary unconditional disturbances by transforming the data appropriately. In this case, take the log of the response series. Difference the logged responses. It is recommended to transform the predictors the same way as the responses to maintain the original interpretation of their relationship. However, do not transform the predictors in this case because they contain negative values. Reestimate the regression model using the transformed responses, and plot the residuals.

dLogY = diff(log(y));
RegMdl2 = fitlm(X(2:end,:),dLogY);

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

h = adftest(RegMdl2.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.