www.gusucode.com > econ 案例源码程序 matlab代码 > econ/ParametricTrendEstimationExample.m
%% Parametric Trend Estimation % This example shows how to estimate nonseasonal and seasonal trend components % using parametric models. The time series is monthly accidental deaths % in the U.S. from 1973 to 1978 (Brockwell and Davis, 2002). % Copyright 2015 The MathWorks, Inc. %% Step 1: Load the Data % Load the accidental deaths data set. load(fullfile(matlabroot,'examples','econ','Data_Accidental.mat')) y = Data; T = length(y); figure plot(y/1000) h1 = gca; h1.XLim = [0,T]; h1.XTick = 1:12:T; h1.XTickLabel = datestr(dates(1:12:T),10); title 'Monthly Accidental Deaths'; ylabel 'Number of Deaths (in thousands)'; hold on %% % The data shows a potential quadratic trend and a strong seasonal component % with periodicity 12. %% Step 2: Fit Quadratic Trend % Fit the polynomial % % $$T_t = \beta_0 + \beta_1t + \beta_2t^2$$ % % to the observed series. t = (1:T)'; X = [ones(T,1) t t.^2]; b = X\y; tH = X*b; h2 = plot(tH/1000,'r','LineWidth',2); legend(h2,'Quadratic Trend Estimate') hold off %% Step 3. Detrend Original Series. % Subtract the fitted quadratic line from the original data. xt = y - tH; %% Step 4. Estimate Seasonal Indicator Variables % Create indicator (dummy) variables for each month. The first indicator % is equal to one for January observations, and zero otherwise. The second % indicator is equal to one for February observations, and zero otherwise. % A total of 12 indicator variables are created for the 12 months. Regress % the detrended series against the seasonal indicators. mo = repmat((1:12)',6,1); sX = dummyvar(mo); bS = sX\xt; st = sX*bS; figure plot(st/1000) title 'Parametric Estimate of Seasonal Component (Indicators)'; h3 = gca; h3.XLim = [0,T]; ylabel 'Number of Deaths (in thousands)'; h3.XTick = 1:12:T; h3.XTickLabel = datestr(dates(1:12:T),10); %% % In this regression, all 12 seasonal indicators are included in the design % matrix. To prevent collinearity, an intercept term is not included (alternatively, % you can include 11 indicators and an intercept term). %% Step 5. Deseasonalize Original Series % Subtract the estimated seasonal component from the original series. dt = y - st; figure plot(dt/1000) title 'Monthly Accidental Deaths (Deseasonalized)'; h4 = gca; h4.XLim = [0,T]; ylabel 'Number of Deaths (in thousands)'; h4.XTick = 1:12:T; h4.XTickLabel = datestr(dates(1:12:T),10); %% % The quadratic trend is much clearer with the seasonal component removed. %% Step 6. Estimate Irregular Component % Subtract the trend and seasonal estimates from the original series. The % remainder is an estimate of the irregular component. bt = y - tH - st; figure plot(bt/1000) title('Irregular Component') h5 = gca; h5.XLim = [0,T]; ylabel 'Number of Deaths (in thousands)'; h5.XTick = 1:12:T; h5.XTickLabel = datestr(dates(1:12:T),10); %% % You can optionally model the irregular component using a stochastic process % model. %% % 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.