2016-11-06 4 views
0

J'ai créé ce code ici:Matlab Kalman Filter Code - Pourquoi ne pas travailler?

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 

charge hon.txt et charge dji.txt

chargent correctement comme 101-by-1 double.

Quand je cours, je reçois le texte suivant:

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


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

Ce qui me surprend parce que dans l'espace de travail, je peux voir que

F = 0 

....

Je suis confus. ..

Répondre

0

Vous devez dire explicitement à MATLAB que lorsque vous dites

resultparam = fminsearch(MyLikelihoodFn, param); 

ce que vous voulez passer dans fminsearch est la fonction MyLikelihoodFn elle-même plutôt que le résultat de l'appeler. (Matlab traite un nom de fonction nue comme une invocation de cette fonction.)

donc mettre un @ devant lui:

resultparam = fminsearch(@MyLikelihoodFn, param); 
+0

fonctionne parfaitement! – user3165675