clear all;
clc;
% Define the global constants
%
% n = number of measured elements (rows of y)
% k = number of states
% s = number of shocks affecting the state
% obs = number of observations (including predicted values)
global g_y g_n g_k g_s g_obs
global g_a g_a0 g_ap g_as g_P g_P0 g_Pp g_LL
global g_Z g_d g_H g_T g_c g_R g_Q
global g_infinity
n=2; g_n=n;
k=2; g_k=k;
s=2; g_s=s;
obs=1000; g_obs=obs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Kalman Filtering
% This is a Multivariate Kalman filtering routine that
% handles missing observations. It is a linear filter, and
% it has been augmented with a fixed interval smoother.
% It is based on that described in Harvey (1989) "Forecasting
% Structural Time Series Models and the Kalman Filter". It
% generally complies with the notation established in that
% text.
% Initialise the global system matrices and filtration
% output matrices
g_y = zeros(n,obs);
g_Z = zeros(k*n,obs);
g_d = zeros(n,obs);
g_H = zeros(n*n,obs);
g_T = zeros(k*k,obs);
g_c = zeros(k,obs);
g_R = zeros(k*s,obs);
g_Q = zeros(s*s,obs);
g_a = zeros(k,obs);
g_a0 = zeros(k,1);
g_ap = a;
g_as = a;
g_P = zeros(k*k,obs);
g_P0 = zeros(k*k,1);
g_Pp = P;
g_Ps = P;
g_LL = zeros(1,obs);
g_infinity = 1e5;
% Define the filtration
%
% The process for each observation except the first:
% (Note that the first observation is the prior!!!)
%
% 1. Form the system matrices for the current period
%
% 2. Get the state and precision matrix predictions for
% the current observation, based on date up to the
% previous observation.
%
% 3. Compute the prediction errors and the prediction
% error variances.
%
% 4. Compute the updated state and precision matrices
% and store them. Also store the prediction error
% variances and the prediction errors themselves.
%
% Note that all of the filtration output is stored in global
% matrices which makes the filter calls faster at the cost
% of some robustness of the code! This means that there
% are no return parameters to worry about.
% Specify the initial values of the hyper parameters
parms = [0.25 % sigmah1
0.25 % sigmah2
0.15 % sigmaq1
0.15] ; % sigmaq2
% Specify the data
% For this example, the data is simulated.
data = zeros(n,obs);
for i=2:obs
a(1,i)=a(1,i-1)+parms(3)*randn(1,1);
a(2,i)=a(2,i-1)+parms(4)*randn(1,1);
y(1,i)=a(1,i-1)+parms(1)*randn(1,1);
y(2,i)=a(2,i-1)+parms(2)*randn(1,1);
end
data=y;
% Define the global matrices used for filtering
% Note that the global variables start with _
%
% y = z * a + d + e
% a = T * a(-1) + c + R * u
%
% Var(e) = H
% Var(u) = Q
%
% Store them vectorised and stacked horizontally.
%
% Thus we have:
%
% For the measurement equation:
% y ~ n*1 by obs
% Z ~ k*n by obs
% d ~ n*1 by obs
% H ~ n*n by obs
%
% For the transition equation:
% T ~ k*k by obs
% c ~ k*1 by obs
% R ~ k*s by obs
% Q ~ s*s by obs
%
% For the state vector we have:
% a ~ k*1 by obs conditioned on the current period
% ap ~ k*1 by obs conditioned on the previous period
% as ~ k*1 by obs conditioned on all of the data
% a0 ~ k by 1 initial state vector
%
% For the state vector precision matrix we have:
% P ~ k*k by obs conditioned on the current period
% Pp ~ k*k by obs conditioned on the previous period
% Ps ~ k*k by obs conditioned on all of the data
% P0 ~ k by k initial state vector precision matrix
%
% Finally we have to store the loglikelihood function
% LL ~ 1 by obs
%
% !!!!!!!!!!!!!!!!!!!!!NOTE!!!!!!!!!!!!!!!!!!!!!!!!
%
% This is where the definition of the model must be
% developed. It needs to be altered to suit each model
% specification.
%
% The example below specifies a model that is a mixture of
% two independent first order autoregressive processes. The
% measurement equation shocks are zero so the mixed AR process
% is not observed with noise. This last point translates into
% having the H matrix just being a set of zeros.
%
% !!!!!!!!!!!!!!!!!!!!!NOTE!!!!!!!!!!!!!!!!!!!!!!!!
%
% The matrices are initialised where possible to
% the values that will be used in the filtering process.
%
% Where the matrices depend on hyper parameters, however,
% these will need to be set just before the filtering
% occurs so that the most recently obtained parameter
% estimates are used by the filter. This is especially
% relevant when engaging in MLE using the Kalman Filter.
y=data;
Z(1,:)=ones(1,obs);
Z(4,:)=ones(1,obs);
R(1,:)=ones(1,obs);
R(4,:)=ones(1,obs);
T(1,:)=ones(1,obs);
T(4,:)=ones(1,obs);
P(1,1)=infinity;
P(4,1)=infinity;
% Call SysMat
[H,Q,a0,P0]=SysMat(parms);
[LL]=KalmanFilter(Z,H,T,R,Q);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%