SNopp

mail@pastecode.io avatar
unknown
matlab
2 years ago
2.3 kB
1
Indexable
Never
%% comp3 2.2 Kalman filtering of time series
addpath '/Users/noahakesson/MATLAB-Drive/Comp2/Tidserie/functions'
addpath '/Users/noahakesson/MATLAB-Drive/Comp2/Tidserie/data'
close all
clear all

%%
% Simulate an AR(2) process.
rng(0)                                          % Set the seed (just done for the lecture!)
N  = 500;
A0_start = -0.95;
A0_stop  = -0.45;
A1 = linspace( A0_start,A0_start, N);
A2 = linspace( A0_start,A0_start, N)
%A1(N/2:end) = A0_stop;                         % Abruptly change a_1.
y = zeros(N,1);
e = randn( N, 1 );
for k=3:N                                       % Implement filter by hand to allow a1 to change.
    y(k) = e(k) - A1(k)*y(k-1) - A2(k)*y(k-2);
end

rng(0)                                          % Set the seed (just done for the lecture!)
N  = 500;
A0_start = -0.95;
A0_stop  = -0.45;
A1 = linspace( A0_start,A0_start, N);
A2 = linspace( A0_start,A0_start, N)
%A1(N/2:end) = A0_stop;                         % Abruptly change a_1.
y = zeros(N,1);
e = randn( N, 1 );
for k=3:N                                       % Implement filter by hand to allow a1 to change.
    y(k) = e(k) - A1(k)*y(k-1) - A2(k)*y(k-2);
end
% Define the state space equations .
A = [ 1 0 ; 0 1 ] ;
Re = [ 0.004 0 ; 0 0 ] ;            % State covariance matrix
Rw = 1.25                           % Observation variance
% Set some initialvalues
Rxx_1 = 10*eye ( 2 ) ;              % Initialstate variance
xtt_1 = [ 0 0 ]';                   % Initialstatevalues
% Vectors to store values in
Xsave = zeros( 2 ,N) ;              % Stored states
ehat = zeros( 1 ,N) ;               % Predict on r e s i d u a l
yt1 = zeros( 1 ,N) ;                % One s t e p p r e d i c t i o n
yt2 = zeros( 1 ,N) ;                % Two s t e p p r e d i c t i o n

% The f i l t e r use data up to t ime t􀀀1 to p r e d i c t v a lue at t ,
% then update using the prediction error. Why do we start
% from t=3? Why stop at N-2?
for t=3:N-2
Ct = [  -y(t-1) -y(t-2)  ]; 
yhat( t ) = ?
ehat( t ) = y(t) - yhat(t-1) 

% Update
Ryy = ? 
Kt = ? 
xtt = ? 
Rxx = ? 
% Pr e d i c t the ne x t s t a t e
xtt_1 = ? 
Rxx_1 = ? 

% Form 2􀀀s t e p p r e d i c t i o n . Ignore t h i s par t at f i r s t .

Ct1 = [ ? ? ]; 
yt1(t+1) = ? 
Ct2 = [ ? ? ] ; 
yt2(t+2) = ? 

% Stor e the s t a t e v e c t o r
Xsave(:, t) = xtt;
end