%% Transient Effects Are Randomly Spread
% This example examines regression lines of regression models with ARMA errors when the transient effects are randomly spread with respect to the joint distribution of the predictor and response.
% Specify the regression model with ARMA(2,1) errors:
% y_t=3 + 2X_t + u_t
% u_t=0.9u_{t-2} - 0.4u_{t-2}+\varepsilon_t + 0.8\varepsilon_{t-1},
% where \varepsilon_t is Gaussian with mean 0 and variance 1. Plot the impulse response function.

Mdl0 = regARIMA('AR',{0.9,-0.4},'MA',{0.8},'Beta',2,...
    'Variance',1,'Intercept',3);

figure
impulse(Mdl0)

% The unconditional disturbances seem to settle after the 10th lag. Therefore, the transient effects end at the 10th lag.
%
% Simulate a univariate, Gaussian predictor series with mean 0 and variance 1. Simulate 100 paths from Mdl0.

rng(5);         % For reproducibility
T = 50;         % Sample size
numPaths = 100; % Number of paths

X = randn(T,1);                                % Full predictor series
Y = simulate(Mdl0,T,'numPaths',numPaths,'X',X); % Full response series

endTrans = 10;
truncX = X((endTrans+1):end);   % Predictor without transient effects    
truncY = Y((endTrans+1):end,:); % Response without transient effects

%
% Fit the model to each simulated response path separately for the full and truncated series.

Mdl = regARIMA(2,0,1); % Empty model for estimation
beta1 = zeros(2,numPaths);
beta2 = beta1;

for i = 1:numPaths
    EstMdl1 = estimate(Mdl,Y(:,i),'X',X,'display','off');
    EstMdl2 = estimate(Mdl,truncY(:,i),'X',truncX,'display','off');
    beta1(:,i) = [EstMdl1.Intercept; EstMdl1.Beta];
    beta2(:,i) = [EstMdl2.Intercept; EstMdl2.Beta];
end

% beta1 is a 2-by-numPaths matrix containing the estimated intercepts and slopes for each simulated data set.  beta2 is a 2-by-numPaths matrix containing the estimated intercepts and slopes for the truncated, simulated data sets.
%
% Compare the simulated regression lines between the full and truncated series. For one of the paths, plot the simulated data and its corresponding regression lines.

betaBar1 = mean(beta1,2);
betaBar2 = mean(beta2,2);
displayresults(betaBar1,betaBar2) % See Supporting Functions section
figure
plot(X,Y(:,1),'.')
hold on
plot(X(1:endTrans),Y(1:endTrans),'ro')
plot([min(X) max(X)],beta1(1,1) + beta1(2,1)*[min(X) max(X)],'b')
plot([min(truncX) max(truncX)],...
    beta2(1,1) + beta2(2,1)*[min(truncX) max(truncX)],'r')
legend('Data','Truncated Data','With Transient Effects',...
    'Without Transient Effects','Location','NorthWest')
xlabel('x')
ylabel('y')
text(0,-3,sprintf('\\beta_0 = %0.4g',beta1(1,1)),'Color',[0,0,1])
text(0,-4,sprintf('\\beta_0 = %0.4g',beta2(1,1)),'Color',[1,0,0])
text(2,-3,sprintf('\\beta_1 = %0.4g',beta1(2,1)),'Color',[0,0,1])
text(2,-4,sprintf('\\beta_1 = %0.4g',beta2(2,1)),'Color',[1,0,0])
hold off

% The table in the Command Window displays the simulation averages of the intercept and slope of the regression model. The results suggest the regression line corresponding to the analysis including the full data set is parallel to the regression line corresponding to the truncated data set. In other words, the slope is mostly unaffected by accounting for transient effects, but the intercept is slightly affected.
%
%% Supporting Functions

function displayresults(b1,b2)
fprintf('Transient Effects  | Sim. Mean of Intercept    | Sim. Mean of Slope\n')
fprintf('===================================================================\n')
fprintf('Include            | %0.6g                   | %0.6g\n',b1(1),b1(2))
fprintf('Without            | %0.6g                   | %0.6g\n',b2(1),b2(2))
end

% Copyright 2015 The MathWorks, Inc.