Search code examples
matlabkalman-filter

Matlab Kalman Filter Code - Why not working?


I created this code here:

clear all;

%  State space reprsentation to be forcasted by kalman filter
%   zhi(t+1) = F*zhi(t) + v(t+1)   --> unbobserved varaibles
%   v~N(0,Q)
%   y(t) = A'*x(t) + H'*zhi(t) + w(t)
%   w~N(0,R)

global y;
global x;
global Hvec;
%%----    Enter  Input parameters
load hon.txt %filename for stock prices
load dji.txt %filename for index prices
n=100; %no. of points to consider
offset=1; %use 1 for daily return or 30 for monthly return
%-------------------------------


datapts=1:offset:(n+1)*offset;
dji=dji(datapts);
hon=hon(datapts);

%Hvec=(dji(1:n)-dji(2:n+1))./dji(2:n+1); %index returns process
%y=(hon(1:n)-hon(2:n+1))./hon(2:n+1); %index returns process

Hvec=log(dji(1:n)./dji(2:n+1)); %index returns process
y=log(hon(1:n)./hon(2:n+1)); %stock returns process

Hvec=flipud(Hvec);
y=flipud(y);


x=ones(n,1);




param=zeros(5,1);
F=0.5;
F=-log(1/F-1);
param(1)=F;
param(2)=0.2;
param(3)=1;
param(4)=0.2;
param(5)=0.5;

resultparam=fminsearch(MyLikelihoodFn,param)


F=resultparam(1)
F=1/(1+exp(-F))
Q=resultparam(2)^2
A=resultparam(3)
R=resultparam(4)^2
betai=resultparam(5)

n=size(y,1);
P=Q;
Ezhi=0.01;
Ezhivec=zeros(n,1);
Ezhivec(1)=Ezhi;

for i=2:n
  yt=y(i);
  xt=x(i);
  H=Hvec(i);
  Ezhi = F*Ezhi + F*P*H*inv(H'*P*H+R)*(yt-A'*xt-H'*betai-H'*Ezhi);
  P = F*P*F' - F*P*H*inv(H'*P*H+R)*H'*P*F' + Q;
  Ezhivec(i)=Ezhi;
end
Ezhivec=Ezhivec+betai;
test=[Ezhivec Hvec y];
tmp=1:n;
subplot(3,1,1);
plot(tmp,y,'-');
subplot(3,1,2);
plot(tmp,Hvec,'-');
subplot(3,1,3);
plot(tmp,Ezhivec,'-');

%plot(tmp,Ezhivec,'-',tmp,Hvec,'-',tmp,y,'-');
%plot(tmp,zhi,'-',tmp,Ezhivec,'-');


%% ------------------
%test=[zhi y]

function ret=MyLikelihoodFn(p)
  global y;
  global x;
  global Hvec;
  F=p(1);
  F=1/(1+exp(-F));
  Q=p(2)^2;
  A=p(3);
  R=p(4)^2;
  betai=p(5);
  n=size(y,1);
  P=Q;
  Ezhi=0;
  Ezhivec=zeros(n,1);
  Ezhivec(1)=Ezhi;
  tmpsum=0;
  tmp1=-(n/2)*log(2*pi);
  for i=2:n
    yt=y(i);
    xt=x(i);
    H=Hvec(i);
    Ezhi = F*Ezhi + F*P*H*inv(H'*P*H+R)*(yt-A'*xt-H'*betai-H'*Ezhi);
    P = F*P*F' - F*P*H*inv(H'*P*H+R)*H'*P*F' + Q;
    Ezhivec(i)=Ezhi;
    tmpmat = H'*P*H + R;
    tmp2 = -0.5*log(det(tmpmat));
    tmpmat2 = yt - A'*xt - H'*betai - H'*Ezhi;
    tmp3=-0.5*tmpmat2'*inv(tmpmat)*tmpmat2;
    tmpsum=tmp1+tmp2+tmp3;
  end
  ret=-tmpsum;   
end

load hon.txt and load dji.txt

are correctly loading as 101-by-1 doubles.

When I run, I get the following:

Error in kalmanbeta>MyLikelihoodFn (line 95)
  F=p(1);


Error in kalmanbeta (line 50)
resultparam=fminsearch(MyLikelihoodFn,param)

Which surprise me because in the workspace I can see that

F = 0

....

I am confused...


Solution

  • You need to tell MATLAB explicitly that when you say

    resultparam = fminsearch(MyLikelihoodFn, param);
    

    what you want to pass into fminsearch is the function MyLikelihoodFn itself rather than the result of calling it. (MATLAB treats a bare function name as an invocation of that function.)

    So put an @ in front of it:

    resultparam = fminsearch(@MyLikelihoodFn, param);