clear all
load data.m;%data starts from december88
pi=data(:,1);%actual inflation
pia=data(:,2);%expected inflation
targetdif=(data(:,3)-pi)^2; %output gap

%a=[0.5 0.25 0.125 0.0625 0.0625];%weights for calculating ro and pip(past inflations) 
%in the loop below we calculate ro(t)and pip(t) starting from december89


%the calculation starts from t0=?
%therefore first calculation of the state vector takes t+1 as ?
dif = pi-pia;
F0 = [0];
sig=0.1;
Q=sig;
R=0.1;
sayHat=zeros(1,60);
sayHat(1,1)=1;
P1giv0=Q;
% vecP1giv0=inv(eye(9)-kron(F0,F0))*vec(Q)
% P1giv0=reshape(vecP1giv0,3,3);
% for t=1:100  %%this loop may be used as a fine tuning mechanism for the Kalman filter...
for i=1:59; %Kalman Filter equations start here....
    F = [1];
    %F = [0 dif(i) dif(i)*rot(i);0 1 rot(i);0 0 1];
    H=[dif];
    if i==1
        sayHat1giv1=sayHat(:,1)+P1giv0*H*inv(H'*P1giv0*H+R)*(pia(2)-H(1)*sayHat(:,1)-pia(1));
        P1giv1=P1giv0-P1giv0*H*inv(H'*P1giv0*H+R)*H'*P1giv0;
        sayHat(:,2)=F*sayHat1giv1;
        PtArti1Givt=F*P1giv1*F'+Q;
    else
        sayHattGivt=sayHat(:,i)+PtArti1Givt*H(i)*inv(H(i)*PtArti1Givt*H(i)+R)*(pi(i+1)-H(i)*sayHat(:,i)-pia(i));
        PtGivt=PtArti1Givt-(PtArti1Givt*H(i))*inv(H(i)*PtArti1Givt*H(i)+R)*H(i)*PtArti1Givt;
        sayHat(:,i+1)=F*sayHattGivt;
        PtArti1Givt=F*PtGivt*F'+Q;
    end
end
%end
plot(sayHat(1,25:60),'b'); %hold on; 
%plot(pi(29:64),'.-g');hold on;plot(pia(29:64),':m');
legend('theta-Kalman');%,'actual pi(t+1)','expected pi(t+1)');
xlabel('# of period');
ylabel('theta');
cor1=corrcoef(sayhat(29:64),targetdif(29:64))
