www.gusucode.com > econ 案例源码程序 matlab代码 > econ/SmoothTimeVaryingStateSpaceModelExample.m
%% Smooth Time-Varying State-Space Model % This example shows how to generate data from a known model, fit a % state-space model to the data, and then smooth the states. %% % Suppose that a latent process comprises an AR(2) and an MA(1) model. % There are 50 periods, and the MA(1) process drops out of the model for % the final 25 periods. Subsequently, the state equation for the first 25 % periods is % % $$\begin{array}{l} % {x_{1,t}} = 0.7{x_{1,t - 1}} - 0.2{x_{1,t - 2}} + {u_{1,t}}\\ % {x_{2,t}} = {u_{2,t}} + 0.6{u_{2,t - 1}}, % \end{array}$$ % % and for the last 25 periods, it is % % $${x_{1,t}} = 0.7{x_{1,t - 1}} - 0.2{x_{1,t - 2}} + {u_{1,t}},$$ % % where $u_{1,t}$ and $u_{2,t}$ are Gaussian with mean 0 and standard % deviation 1. %% % Assuming that the series starts at 1.5 and 1, respectively, generate a % random series of 50 observations from $x_{1,t}$ and $x_{2,t}$. % Copyright 2015 The MathWorks, Inc. T = 50; ARMdl = arima('AR',{0.7,-0.2},'Constant',0,'Variance',1); MAMdl = arima('MA',0.6,'Constant',0,'Variance',1); x0 = [1.5 1; 1.5 1]; rng(1); x = [simulate(ARMdl,T,'Y0',x0(:,1)),... [simulate(MAMdl,T/2,'Y0',x0(:,2));nan(T/2,1)]]; %% % The last 25 values for the simulated MA(1) data are |NaN| values. %% % Suppose further that the latent processes are measured using % % $${y_t} = 2\left( {{x_{1,t}} + {x_{2,t}}} \right) + {\varepsilon _t},$$ % % for the first 25 periods, and % % $${y_t} = 2{x_{1,t}} + {\varepsilon _t}$$ % % for the last 25 periods, where $\varepsilon_t$ is Gaussian with mean 0 % and standard deviation 1. % % Use the random latent state process (|x|) and the observation equation to % generate observations. y = 2*nansum(x')'+randn(T,1); %% % Together, the latent process and observation equations compose a % state-space model. Supposing that the coefficients are unknown % parameters, the state-space model is % % $$\begin{array}{l} % \left[ {\begin{array}{*{20}{c}} % {{x_{1,t}}}\\ % {{x_{2,t}}}\\ % {{x_{3,t}}}\\ % {{x_{4,t}}} % \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} % {{\phi _1}}&{{\phi _2}}&0&0\\ % 1&0&0&0\\ % 0&0&0&{{\theta _1}}\\ % 0&0&0&0 % \end{array}} \right]\left[ {\begin{array}{*{20}{c}} % {{x_{1,t - 1}}}\\ % {{x_{2,t - 1}}}\\ % {{x_{3,t - 1}}}\\ % {{x_{4,t - 1}}} % \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} % 1&0\\ % 0&0\\ % 0&1\\ % 0&1 % \end{array}} \right]\left[ {\begin{array}{*{20}{c}} % {{u_{1,t}}}\\ % {{u_{2,t}}} % \end{array}} \right]\\ % {y_t} = a({x_{1,t}} + {x_{3,t}}) + {\varepsilon _t} % \end{array}$$ % % for the first 25 periods, % % $$\begin{array}{l} % \left[ {\begin{array}{*{20}{c}} % {{x_{1,t}}}\\ % {{x_{2,t}}} % \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} % {{\phi _1}}&{{\phi _2}}&0&0\\ % 1&0&0&0 % \end{array}} \right]\left[ {\begin{array}{*{20}{c}} % {{x_{1,t - 1}}}\\ % {{x_{2,t - 1}}}\\ % {{x_{3,t - 1}}}\\ % {{x_{4,t - 1}}} % \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} % 1\\ % 0 % \end{array}} \right]{u_{1,t}}\\ % {y_t} = b{x_{1,t}} + {\varepsilon _t} % \end{array}$$ % % for period 26, and % % $$\begin{array}{l} % \left[ {\begin{array}{*{20}{c}} % {{x_{1,t}}}\\ % {{x_{2,t}}} % \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} % {{\phi _1}}&{{\phi _2}}\\ % 1&0 % \end{array}} \right]\left[ {\begin{array}{*{20}{c}} % {{x_{1,t - 1}}}\\ % {{x_{2,t - 1}}} % \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} % 1\\ % 0 % \end{array}} \right]{u_{1,t}}\\ % {y_t} = b{x_{1,t}} + {\varepsilon _t} % \end{array}$$ % % for the last 24 periods. %% % Write a function that specifies how the parameters in |params| map to the % state-space model matrices, the initial state values, and the type of % state. % % <include>AR2MAParamMap.m</include> % %% % Save this code as a file named |AR2MAParamMap| on your MATLAB(R) path. %% % Create the state-space model by passing the function |AR2MAParamMap| as a % function handle to |ssm|. Mdl = ssm(@(params)AR2MAParamMap(params,T)); %% % |ssm| implicitly creates the state-space model. Usually, you % cannot verify an implicitly defined state-space model. %% % Pass the observed responses (|y|) to |estimate| to estimate the % parameters. Specify an arbitrary set of positive initial values for the % unknown parameters. params0 = 0.1*ones(5,1); EstMdl = estimate(Mdl,y,params0); %% % |EstMdl| is an |ssm| model containing the estimated coefficients. % Likelihood surfaces of state-space models might contain local maxima. % Therefore, it is good practice to try several initial parameter values, % or consider using |refine|. %% % Smooth the states and estimate the variance-covariance matrices of the % smoothed states by passing |EstMdl| and the observed responses to % |smooth|. [~,~,Output]= smooth(EstMdl,y); %% % |Output| is a |T|-by-1 structure array containing the smoothed states and % their variance-covariance matrices, among other things. %% % Extract the smoothed states and their variances from the cell arrays. % Recall that the two, different states are in positions 1 and 3. The % states in positions 2 and 4 help specify the processes of interest. stateIndx = [1 3]; % State Indices of interest SmoothedStates = NaN(T,numel(stateIndx)); SmoothedStatesCov = NaN(T,numel(stateIndx)); for t = 1:T maxInd1 = size(Output(t).SmoothedStates,1); maxInd2 = size(Output(t).SmoothedStatesCov,1); mask1 = stateIndx <= maxInd1; mask2 = stateIndx <= maxInd2; SmoothedStates(t,mask1) = ... Output(t).SmoothedStates(stateIndx(mask1),1); SmoothedStatesCov(t,mask2) = ... diag(Output(t).SmoothedStatesCov(stateIndx(mask2),... stateIndx(mask2))); end %% % Plot the true state values, the smoothed state values, and their % individual 95% Wald-type confidence intervals for each model. AR2SSCIlb = SmoothedStates(:,1) - 1.95*sqrt(SmoothedStatesCov(:,1)); AR2SSCIub = SmoothedStates(:,1) + 1.95*sqrt(SmoothedStatesCov(:,1)); AR2SSIntervals = [AR2SSCIlb AR2SSCIub]; MA1SSCIlb = SmoothedStates(:,2) - 1.95*sqrt(SmoothedStatesCov(:,2)); MA1SSCIub = SmoothedStates(:,2) + 1.95*sqrt(SmoothedStatesCov(:,2)); MA1SSIntervals = [MA1SSCIlb MA1SSCIub]; figure plot(1:T,x(:,1),'-k',1:T,SmoothedStates(:,1),':r',... 1:T,AR2SSIntervals,'--b','LineWidth',2); title('AR(2) State Values') xlabel('Period') ylabel('State Value') legend({'True state values','Smoothed state values',... '95% Confidence Intervals'}); figure plot(1:T,x(:,2),'-k',1:T,SmoothedStates(:,2),':r',... 1:T,MA1SSIntervals,'--b','LineWidth',2); title('MA(1) State Values') xlabel('Period') ylabel('State Value') legend({'True state values','Smoothed state values',... '95% Confidence Intervals'});