www.gusucode.com > econ 案例源码程序 matlab代码 > econ/EstimateaRegressionModelwithMultiplicativeARIMAErrorsExample.m

    %% Estimate a Regression Model with Multiplicative ARIMA Errors  
% This example shows how to fit a regression model with multiplicative ARIMA
% errors to data using |estimate|.   

% Copyright 2015 The MathWorks, Inc.


%% 
% Load the Airline data set from the MATLAB(R) root folder, and load the
% Recession data set. Plot the monthly passenger totals and the log of the
% totals. 
load(fullfile(matlabroot,'examples','econ','Data_Airline.mat'))
load Data_Recessions

y = Data;
logY = log(y);

figure
subplot(2,1,1)
plot(y)
title('{\bf Monthly Passenger Totals (Jan1949 - Dec1960)}')
datetick
subplot(2,1,2)
plot(log(y))
title('{\bf Monthly Passenger Log-Totals (Jan1949 - Dec1960)}')
datetick    

%%
% The log transformation seems to linearize the time series.  

%% 
% Construct the predictor (|X|), which is whether the country was in a recession
% during the sampled period. A 0 in row _t_ means the country was not in
% a recession in month _t_, and a 1 in row _t_ means that it was in a recession
% in month _t_. 
X = zeros(numel(dates),1); % Preallocation
for j = 1:size(Recessions,1)
    X(dates >= Recessions(j,1) & dates <= Recessions(j,2)) = 1;
end  

%% 
% Fit the simple linear regression model
%
% $${y_t} = c + {X_t}\beta  + {u_t}$$
%
% to the data. 
Fit = fitlm(X,logY); 

%%
% |Fit| is a |LinearModel| that contains the least squares estimates.  

%% 
% Check for standard linear model assumption departures by plotting the
% residuals several ways. 
figure
subplot(2,2,1)
plotResiduals(Fit,'caseorder','ResidualType','Standardized',...
    'LineStyle','-','MarkerSize',0.5)
subplot(2,2,2)
plotResiduals(Fit,'lagged','ResidualType','Standardized')
subplot(2,2,3)
plotResiduals(Fit,'probability','ResidualType','Standardized')
subplot(2,2,4)
plotResiduals(Fit,'histogram','ResidualType','Standardized')

r = Fit.Residuals.Standardized;
figure
subplot(2,1,1)
autocorr(r)
subplot(2,1,2)
parcorr(r)    

%%
% The residual plots indicate that the unconditional disturbances are autocorrelated.
% The probability plot and histogram seem to indicate that the unconditional
% disturbances are Gaussian.    

%%
% The ACF of the residuals confirms that the unconditional disturbances
% are autocorrelated.  

%% 
% Take the 1st difference of the residuals and plot the ACF and PACF of
% the differenced residuals. 
dR = diff(r);

figure
subplot(2,1,1)
autocorr(dR,50)
subplot(2,1,2)
parcorr(dR,50)    

%%
% The ACF shows that there are significantly large autocorrelations, particularly
% at every 12th lag. This indicates that the unconditional disturbances
% have 12th degree seasonal integration.  

%% 
% Take the first and 12th differences of the residuals. Plot the differenced
% residuals, and their ACF and PACF. 
DiffPoly = LagOp([1 -1]);
SDiffPoly = LagOp([1 -1],'Lags',[0, 12]);
diffR = filter(DiffPoly*SDiffPoly,r);

figure
subplot(2,1,1)
plot(diffR)
axis tight
subplot(2,2,3)
autocorr(diffR)
axis tight
subplot(2,2,4)
parcorr(diffR)
axis tight    

%%
% The residuals resemble white noise (with possible heteroscedasticity).
% According to Box and Jenkins (1994), Chapter 9, the ACF and PACF indicate
% that the unconditional disturbances are an
% $\rm{IMA}(0,1,1)\times(0,1,1)_{12}$ model.

%% 
% Specify the regression model with $\rm{IMA}(0,1,1)\times(0,1,1)_{12}$
% errors: 
%
% $$\begin{array}{c}{y_t} = {X_t}\beta  + {u_t}\\\left( {1 - L} \right)\left( {1 - {L^{12}}} \right){u_t} = \left( {1 + {b_1}L} \right)\left( {1 + {B_{12}}{L^{12}}} \right){\varepsilon _t}.\end{array}$$
%
Mdl = regARIMA('MALags',1,'D',1,'Seasonality',12,'SMALags',12)  

%% 
% Partition the data set into the presample and estimation sample so that
% you can initialize the series. |P| = |Q| = 13, so the presample should
% be at least 13 periods long. 
preLogY = logY(1:13);   % Presample responses
estLogY = logY(14:end); % Estimation sample responses
preX = X(1:13);         % Presample predictors
estX = X(14:end);       % Estimation sample predictors  

%% 
% Obtain presample unconditional disturbances from a linear regression of
% the presample data. 
PreFit = fitlm(preX,preLogY);...
    % Presample fit for presample residuals
EstFit = fitlm(estX,estLogY);...
    % Estimation sample fit for the intercept
U0 = PreFit.Residuals.Raw;  

%% 
% If the error model is integrated, then the regression model intercept
% is not identifiable. Set |Intercept| to the estimated intercept from a
% linear regression of the estimation sample data. Estimate the regression
% model with IMA errors. 
Mdl.Intercept = EstFit.Coefficients{1,1};
EstMdl = estimate(Mdl,estLogY,'X',estX,'U0',U0);  

%% 
% |MA{1}| and |Beta1| are not significantly different from 0. You can remove
% these parameters from the model, possibly add other parameters (e.g.,
% AR parameters), and compare multiple model fits using |aicbic|. Note that
% the estimation and presample should be the same over competing models.
%%
% References:
%
% Box, G. E. P., G. M. Jenkins, and G. C. Reinsel. _Time Series Analysis: Forecasting and Control_. 3rd ed. Englewood Cliffs, NJ: Prentice Hall, 1994.