close all
clear all
clc

cz_data;

% cas
% t ... roky
[n,p] = size(data);
tts = t(1)+4/4;    % start. rok
ttf = t(end);      % konec. rok
tt  = tts:1/4:ttf; % orezany vektor let
ts  = 5;           % time start
n   = length(tt);  % pocet dat

% MODELOVA DATA
cpi     = CPI(ts:end);
cpi4    = CPI4(ts:end);
cpi4_4  = CPI4(1:end-ts+1);
rs4     = RS4(ts:end);
gr_rs4  = GR_RS4(ts:end);
koruna_dm = SS_DM(ts:end);
gr_cpi    = GR_CPI(ts:end);
gr_cpi_4  = GR_CPI(1:end-ts+1);
lz_eq1 = LZ_EQ(ts:end);
prem   = PREM(ts:end);
% Priprava dat - vypoctete

% inflace
pie = 100*(log(cpi) - log(cpi4_4));

% inflace zahranicni
gr_pie = 100*(log(gr_cpi) - log(gr_cpi_4));

% realne sazby
rr4 = rs4 - pie;
gr_rr4 = gr_rs4 - gr_pie;

% realny kurz
ls_dm = 100*log(1./koruna_dm);
lz    = ls_dm + 100*log(cpi) - 100*log(gr_cpi);

%% Vypoctete skoln LZ, neboli prolozte LZ primku jejiz parametry odhadnete pomoci OLS
b_ = [ ones(size(lz)), (1:length(lz))' ] \ lz;
DRIFT = 4 * b_(2);

% model

%% Measurement equations

%% rr4(t)    = rr4_eq(t)          + w1(t)
%% gr_rr4(t) = gr_rr4_eq(t)       + w2(t)
%% lz(t)     = lz_eq(t)           + w3(t)
%% prem      = prem_eq            + w4(t)

%% State transition
%% rr4_eq(t)      = rr4_eq(t-1) + v1(t)
%% gr_rr4_eq(t)   = gr_rr4_eq(t-1) + v2(t)
%% lz_eq(t)       = lz_eq(t-1)+ DRIFT/4  + v3(t)
%% 4*lz_eq(t)     = 4*lz_eq(t+1) + rr4_eq(t) - gr_rr4_eq(t) - prem_eq(t) + v4(t) ==> tady budeme pocitat prem_eq

%% Matrix denotation

%% Y = C*X + D*U + w
%% X = A*X(-1) + B*U + v


% Prepis modelu do stavoveho tvaru

%A0
%      rr4_eq(t) gr_rr4_eq  lz_eq(t)   prem_eq
A0 = [ 0          0            0             0; %rr4_eq
       0          0            0             0; % gr_rr4_eq
       0          0            0             0; % lz_eq
       1         -1           -4             0]; % prem_eq
   
% A1
%      rr4_eq(t-1) gr_rr4_eq(t-1)  lz_eq(t-1)   prem_eq(t-1)
A1 = [ 1                0            0             0; %rr4_eq
       0                1            0             0; % gr_rr4_eq
       0                0            1             0; % lz_eq
       0                0            0             0]; % prem_eq

% B

B0 = [0 0;
     0 0;
     0 1;
     4 0];
 
 % C
 C = [ 1 0 0 0;
       0 1 0 0;
       0 0 1 0;
       0 0 0 1];
   
 
 % D
 D = [0 0;
      0 0;
      0 0;
      0 0];
  
  % matice
  A = inv(eye(4)-A0)*A1;
  B= inv(eye(4)-A0)*B0;
  
  
  % lz_eq je vysledkem behu HP filtru na lz a prem je spocteno z nominalni verze UIP
  % vystupy
  Y = [ rr4 gr_rr4 lz prem]';
  
  % vstupy
  U = [ lz_eq1 ones(n,1)*DRIFT./4]';  
  
  % pocatecni stavy
  x0 = [rr4(1), gr_rr4(1), lz(1) prem(1)]';
  
  % volte kovariancni matice na zaklade vysledku tak aby odpovidalz vasi ekonomicke predstave,
  % zjistete citlivost KF na kovariancni matice a pocatecni podminky 
  sigma_v = [
        0.1 0 0 0;
        0 0.09 0 0;
        0 0 0.2 0;
        0 0 0 0.9];
  
  sigma_w = [
      5 0 0 0;
      0 5 0 0;
      0 0 5 0;
      0 0 0 5];
  
  sigma_x0 = 0.9*sigma_v;
  
  [mi_p,sigma_p,mi_f,sigma_f,mi_s,sigma_s]=kfs(Y,U,A,B,C,D,sigma_v,sigma_w,x0,sigma_x0);
  
  
  for i=1:n;
      s_rr4_eq(i) = mi_s(1,i);
      p_rr4_eq(i) = mi_p(1,i);
      f_rr4_eq(i) = mi_f(1,i);
      s_gr_rr4_eq(i) = mi_s(2,i);
      p_gr_rr4_eq(i) = mi_p(2,i);
      f_gr_rr4_eq(i) = mi_f(2,i);
      s_lz_eq(i) = mi_s(3,i);
      p_lz_eq(i) = mi_p(3,i);
      f_lz_eq(i) = mi_f(3,i);
      
  end

    
  figure
  plot(tt,rr4,'*-',tt,s_rr4_eq,'.-',tt,f_rr4_eq,':',tt,p_rr4_eq,':')
  legend('data (rr4)','smoothing','filtrace','predikce')
  title('ceska realna urokova mira (rr4) a jeji rovnovazna uroven')
  xlabel('roky')
  ylabel('%')
  
  figure
  plot(tt,gr_rr4,'*-',tt,s_gr_rr4_eq,'.-',tt,f_gr_rr4_eq,':',tt,p_gr_rr4_eq,':')
  legend('data (gr rr4)','smoothing','filtrace','predikce')
  title('nemecka realna urokova mira a jeji rovnovazna uroven')
  xlabel('roky')
  ylabel('%')
  
  figure
  plot(tt,lz./100,'*-',tt,s_lz_eq./100,'.-',tt,f_lz_eq./100,':',tt,p_lz_eq./100,':')
  legend('data (lz/100)','smoothing','filtrace','predikce','location','best')
  title('logaritmus realneho kurzu Kc/DEM a jeho rovnovazna uroven')
  xlabel('roky')
  ylabel('log(Kc/DEM)')