Main Content

Estimate Random Parameter of State-Space Model

This example shows how to estimate a random, autoregressive coefficient of a state in a state-space model. That is, this example takes a Bayesian view of state-space model parameter estimation by using the "zig-zag" estimation method.

Suppose that two states (x1,t and x2,t) represent the net exports of two countries at the end of the year.

  • x1,t is a unit root process with a disturbance variance of σ12.

  • x2,t is an AR(1) process with an unknown, random coefficient and a disturbance variance of σ22.

  • An observation (yt) is the exact sum of the two net exports. That is, the net exports of the individual states are unknown.

Symbolically, the true state-space model is


Simulate Data

Simulate 100 years of net exports from:

  • A unit root process with a mean zero, Gaussian noise series that has variance 0.12.

  • An AR(1) process with an autoregressive coefficient of 0.6 and a mean zero, Gaussian noise series that has variance 0.22.

  • x1,0=x2,0=0.

  • Create an observation series by summing the two net exports per year.

rng(100); % For reproducibility
T = 150;
sigma1 = 0.1;
sigma2 = 0.2;
phi = 0.6;

u1 = randn(T,1)*sigma1;
x1 = cumsum(u1);
Mdl2 = arima('AR',phi,'Variance',sigma2^2,'Constant',0);
x2 = simulate(Mdl2,T,'Y0',0);
y = x1 + x2;

plot([x1 x2 y])
ylabel('Net exports');

Figure contains an axes object. The axes object with xlabel Period, ylabel Net exports contains 3 objects of type line. These objects represent x_1, x_2, y.

The Zig-Zag Estimation Method

Treat ϕ as if it is unknown and random, and use the zig-zag method to recover its distribution. To implement the zig-zag method:

1. Choose an initial value for ϕ in the interval (-1,1), and denote it ϕz.

2. Create the true state-space model, that is, an ssm model object that represents the data-generating process.

3. Use the simulation smoother (simsmooth) to draw a random path from the distribution of the second smoothed states. Symbolically, x2,z,tP(x2,t|yt,x1,t,ϕ=ϕz).

4. Create another state-space model that has this form


In words, ϕz,t is a static state and x2,z,t is an "observed" series with time varying coefficient Ct=x2,z,t-1.

5. Use the simulation smoother to draw a random path from the distribution of the smoothed ϕz,t series. Symbolically, ϕz,tP(ϕ|x2,z,t), where x2,z,t encompasses the structure of the true state-space model and the observations. ϕz,t is static, so you can just reserve one value (ϕz).

6. Repeat steps 2 - 5 many times and store ϕz each iteration.

7. Perform diagnostic checks on the simulation series. That is, construct:

  • Trace plots to determine the burn in period and whether the Markov chain is mixing well.

  • Autocorrelation plots to determine how many draws need removing to obtain a well-mixed Markov chain.

8. The remaining series represents draws from the posterior distribution of ϕ. You can compute descriptive statistics, or plot a histogram to determine the qualities of the distribution.

Estimate Random Coefficient Using Zig-Zag Method

Specify initial values, preallocate, and create the true state-space model.

phi0 = -0.3;             % Initial value of phi
Z = 1000;                % Number of times to iterate the zig-zag method
phiz = [phi0; nan(Z,1)]; % Preallocate

A = [1 0; 0 NaN];
B = [sigma1; sigma2];
C = [1 1];
Mdl = ssm(A,B,C,'StateType',[2; 0]);

Mdl is an ssm model object. The NaN acts as a placeholder for ϕ.

Iterate steps 2 - 5 of the zig-zag method.

for j = 2:(Z + 1);
    % Draw a random path from smoothed x_2 series.
    xz = simsmooth(Mdl,y,'Params',phiz(j-1));
    % The second column of xz is a draw from the posterior distribution of x_2.
    % Create the intermediate state-space model.
    Az = 1;
    Bz = 0;
    Cz = num2cell(xz((1:(end - 1)),2));
    Dz = sigma2;
    Mdlz = ssm(Az,Bz,Cz,Dz,'StateType',2);
    % Draw a random path from the smoothed phiz series.
    phizvec = simsmooth(Mdlz,xz(2:end,2));
    phiz(j) = phizvec(1);
    % phiz(j) is a draw from the posterior distribution of phi

phiz is a Markov chain. Before analyzing the posterior distribution of ϕ, you should assess whether to impose a burn-in period, or the severity of the autocorrelation in the chain.

Determine Quality of Simulation

Draw a trace plot for the first 100, 500, and all of the random draws.

vec = [100 500 Z];
for j = 1:3;
    title('Trace Plot for \phi');
    xlabel('Simulation number');
    axis tight;

Figure contains 3 axes objects. Axes object 1 with title Trace Plot for phi, xlabel Simulation number contains an object of type line. Axes object 2 with title Trace Plot for phi, xlabel Simulation number contains an object of type line. Axes object 3 with title Trace Plot for phi, xlabel Simulation number contains an object of type line.

According to the first plot, transient effects die down after about 20 draws. Therefore, a short burn-in period should suffice. The plot of the entire simulation shows that the series settles around a center.

Plot the autocorrelation function of the series after removing the first 20 draws.

burnOut = 21:Z;


Figure contains an axes object. The axes object with title Sample Autocorrelation Function, xlabel Lag, ylabel Sample Autocorrelation contains 4 objects of type stem, line, constantline.

The autocorrelation function dies out rather quickly. It doesn't seem like autocorrelation in the chain is an issue.

Determine qualities of the posterior distribution of ϕ by computing simulation statistics and by plotting a histogram of the reduced set of random draws.

xbar = mean(phiz(burnOut))
xbar = 0.5104
xstd = std(phiz(burnOut))
xstd = 0.0988
ci = norminv([0.025,0.975],xbar,xstd); % 95% confidence interval

h = gca;
hold on;
simX = linspace(h.XLim(1),h.XLim(2),100);
simPDF = normpdf(simX,xbar,xstd); 
h1 = plot([xbar xbar],h.YLim,'r','LineWidth',2);
h2 = plot([0.6 0.6],h.YLim,'g','LineWidth',2);
h3 = plot([ci(1) ci(1)],h.YLim,'--r',...
    [ci(2) ci(2)],h.YLim,'--r','LineWidth',2);
legend([h1 h2 h3(1)],{'Simulation Mean','True Mean','95% CI'});
h.XTick = sort([h.XTick xbar]);
h.XTickLabel{h.XTick == xbar} = xbar;
h.XTickLabelRotation = 90;

Figure contains an axes object. The axes object contains 6 objects of type histogram, line. These objects represent Simulation Mean, True Mean, 95% CI.

The posterior distribution of ϕ is roughly normal with mean and standard deviation approximately 0.51 and 0.1, respectively. The true mean of ϕ is 0.6, and it is less than one standard deviation to the right of the simulation mean.

Compute the maximum likelihood estimate of ϕ. That is, treat ϕ as a fixed, but unknown parameter, and then estimate Mdl using the Kalman filter and maximum likelihood.

[~,estParams] = estimate(Mdl,y,phi0)
Method: Maximum likelihood (fminunc)
Sample size: 150
Logarithmic  likelihood:     -10.1434
Akaike   info criterion:      22.2868
Bayesian info criterion:      25.2974
      |     Coeff       Std Err     t Stat       Prob  
 c(1) |  0.53590       0.19183    2.79360      0.00521 
      |   Final State   Std Dev    t Stat        Prob  
 x(1) | -0.85059       0.00000   -6.45811e+08   0      
 x(2) |  0.00454        0          Inf          0      
estParams = 0.5359

The MLE of ϕ is 0.54. Both estimates are within one standard deviation or standard error from the true value of ϕ.

See Also

| | | |

Related Examples

More About