www.gusucode.com > econ 案例源码程序 matlab代码 > econ/SmoothTimeVaryingDiffuseStateSpaceModelExample.m
%% Smooth Time-Varying Diffuse State-Space Model % This example shows how to generate data from a known model, fit a % diffuse 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. Consequently, 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. %% % 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 make up a % state-space model. If 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>diffuseAR2MAParamMap.m</include> % %% % Save this code as a file named |diffuseAR2MAParamMap| on your MATLAB(R) % path. %% % Create the diffuse state-space model by passing the function % |diffuseAR2MAParamMap| as a function handle to |dssm|. Mdl = dssm(@(params)diffuseAR2MAParamMap(params,T)); %% % |dssm| implicitly creates the diffuse state-space model. Usually, you % cannot verify diffuse state-space models that are implicitly created. %% % To estimate the parameters, pass the observed responses (|y|) to % |estimate|. 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 a |dssm| model containing the estimated coefficients. % Likelihood surfaces of state-space models might contain local maxima. % Therefore, try several initial parameter values, or consider using % |refine|. %% % Smooth the states and obtain the smoothed state covariance matrix per % period by passing |EstMdl| and the observed responses to |smooth|. [~,~,Output]= smooth(EstMdl,y); %% % |Output| is a |T|-by-1 structure array that contains the smoothed states. %% % Convert |Output| to a table. OutputTbl = struct2table(Output); OutputTbl(1:10,1:4) % Display first ten rows of first four variables %% % The first two rows of the table contain empty cells or zeros, which % correspond to the observations required to initialize the diffuse Kalman % filter. That is, |SwitchTime| is 2. SwitchTime = 2; %% % Extract the smoothed states from |Output|, and compute their 95% % individual, Wald-type confidence intervals. Recall that the two % different states are in positions 1 and 3. The states in positions 2 and % 4 help to specify the processes of interest. stateIdx = [1 3]; % State indices of interest SmoothedStates = nan(T,numel(stateIdx)); CI = nan(T,2,numel(stateIdx)); for t = (SwitchTime + 1):T maxInd = size(Output(t).SmoothedStates,1); mask = stateIdx <= maxInd; SmoothedStates(t,mask) = Output(t).SmoothedStates(stateIdx(mask),1); CovX = Output(t).SmoothedStatesCov(stateIdx(mask),stateIdx(mask)); CI(t,:,1) = SmoothedStates(t,1) + 1.96*sqrt(CovX(1,1))*[-1 1]; if (max(stateIdx(mask)) > 1) CI(t,:,2) = SmoothedStates(t,2) + 1.96*sqrt(CovX(2,2))*[-1 1]; end end SmoothedStates(1:SwitchTime,:) = 0; CI(1:SwitchTime,:,:) = 0; %% % Plot the true state values, the smoothed states, and the 95% % smoothed-state confidence intervals for each model. figure plot(1:T,x(:,1),'b',1:T,SmoothedStates(:,1),'k',1:T,CI(:,:,1),'--r'); title('AR(2) State Values') xlabel('Period') ylabel('State Value') legend({'True state values','Smoothed state values','95% CI'}); figure plot(1:T,x(:,2),'b',1:T,SmoothedStates(:,2),'k',1:T,CI(:,:,2),'--r'); title('MA(1) State Values') xlabel('Period') ylabel('State Value') legend({'True state values','Smoothed state values','95% CI'});