www.gusucode.com > econ 案例源码程序 matlab代码 > econ/SimulateanARMAErrorModelExample.m
%% Simulate an ARMA Error Model % This example shows how to simulate responses from a regression model with % ARMA errors without specifying a presample. % Copyright 2015 The MathWorks, Inc. %% % Specify the regression model with ARMA(2,1) errors: % % $$\begin{array}{l}{y_t} = 2 + {X_t}\left[ \begin{array}{l} - 2\\1.5\end{array} \right] + {u_t}\\{u_t} = 0.9{u_{t - 1}} - 0.1{u_{t - 2}} + {\varepsilon _t} + 0.5{\varepsilon _{t - 1}},\end{array}$$ % % where $\varepsilon_t$ is distributed with 15 degrees of freedom and variance 1. Beta = [-2; 1.5]; Intercept = 2; a1 = 0.9; a2 = -0.1; b1 = 0.5; Variance = 1; Distribution = struct('Name','t','DoF',15); Mdl = regARIMA('AR',{a1, a2},'MA',b1,... 'Distribution',Distribution,'Intercept',Intercept,... 'Beta',Beta,'Variance',Variance); %% % Generate two length _T_ = 50 predictor series by random selection from % the standard Gaussian distribution. T = 50; rng(6); % For reproducibility X = randn(T,2); %% % The software treats the predictors as nonstochastic series. %% % Generate and plot one sample path of responses from |Mdl|. rng(7); ySim = simulate(Mdl,T,'X',X); figure plot(ySim) title('{\bf Simulated Response Series}') %% % |simulate| requires: % % * |P = 2| presample unconditional disturbances to initialize the autoregressive % component of the error series. % * |Q = 1| presample innovations to initialize the moving average component % of the error series. % %% % Without them, as in this case, |simulate| sets the necessary presample % errors to 0. %% % Alternatively, use |filter| to filter a random innovation series through % |Mdl|. rng(7); e = randn(T,1); yFilter = filter(Mdl,e,'X',X); figure plot(yFilter) title('{\bf Simulated Response Series Using Filtered Innovations}') %% % The plots suggest that the simulated responses and the responses generated % from the filtered innovations are equivalent. %% % Simulate 1000 response paths from |Mdl|. Assess transient effects by plotting % the unconditional disturbance (|U|) variances across the simulated paths % at each period. numPaths = 1000; [Y,~,U] = simulate(Mdl,T,'NumPaths',numPaths,'X',X); figure h1 = plot(Y,'Color',[.85,.85,.85]); title('{\bf 1000 Simulated Response Paths}') hold on h2 = plot(1:T,Intercept+X*Beta,'k--','LineWidth',2); legend([h1(1),h2],'Simulated Path','Mean') hold off figure h1 = plot(var(U,0,2),'r','LineWidth',2); hold on theoVarFix = Variance*(a1*b1*(1+a2)+(1-a2)*(1+a1*b1+b1^2))/... ((1+a2)*((1-a2)^2-a1^2)); h2 = plot([1 T],[theoVarFix theoVarFix],'k--','LineWidth',2); title('{\bf Unconditional Disturbance Variance}') legend([h1,h2],'Simulation Variance','Theoretical Variance',... 'Location','Best') hold off %% % The simulated paths follow their theoretical mean, $c + X\beta$, which is % not constant over time (and might look nonstationary). %% % The variance of the process is not constant, but levels off at the theoretical % variance by the 10th period. The theoretical variance of the ARMA(2,1) % error model is: % % $$\frac{{\sigma _\varepsilon ^2\left[ {{a_1}{b_1}\left( {1 + {a_2}} \right) + \left( {1 - {a_2}} \right)\left( {1 + {a_1}{b_1} + b_1^2} \right)} \right]}}{{{{\left( {1 + {a_2}} \right)}^2}\left[ {{{\left( {1 - {a_2}} \right)}^2} - a_1^2} \right]}} = \frac{{\left[ {0.9(0.5)\left( {1 - 0.1} \right) + \left( {1 + 0.1} \right)\left( {1 + 0.9(0.5) + {{0.5}^2}} \right)} \right]}}{{{{\left( {1 - 0.1} \right)}^2}\left[ {{{\left( {1 + 0.1} \right)}^2} - {{0.9}^2}} \right]}} = {\rm{6}}{\rm{.32}}.$$ % %% % You can reduce transient effects by partitioning the simulated data into % a burn-in portion and a portion for analysis. Do not use the burn-in portion % for analysis. Include enough periods in the burn-in portion to overcome % the transient effects. burnIn = 1:10; notBurnIn = burnIn(end)+1:T; Y = Y(notBurnIn,:); X = X(notBurnIn,:); U = U(notBurnIn,:); figure h1 = plot(notBurnIn,Y,'Color',[.85,.85,.85]); hold on h2 = plot(notBurnIn,Intercept+X*Beta,'k--','LineWidth',2); title('{\bf 1000 Simulated Response Paths for Analysis}') legend([h1(1),h2],'Simulated Path','Mean') axis tight hold off figure h1 = plot(notBurnIn,var(U,0,2),'r','LineWidth',2); hold on h2 = plot([notBurnIn(1) notBurnIn(end)],... [theoVarFix theoVarFix],'k--','LineWidth',2); title('{\bf Converged Unconditional Disturbance Variance}') legend([h1,h2],'Simulation Variance','Theoretical Variance') axis tight hold off %% % Unconditional disturbance simulation variances fluctuate around the theoretical % variance due to Monte Carlo sampling error. Be aware that the exclusion % of the burn-in sample from analysis reduces the effective sample size.