%% run the baseline model with the alternative data %%
 
clear all;
load Data;

% use the appropriate data: sample period, (Feb 1972 - June 2017)
U = U_hs;  % Hornstein (2012)'s alternative data
Y = U(48:end,1:5)/1000; % Jan1976 - Jun2017

Si = size(Y);n=Si(:,2);T=Si(:,1);nt=47;

global T;global Y;global U;global nt;
x0 = [0.1;1;-2;0.1;0.01];

fun1 = @ss_init1;
[xmin1,fmin1,flag1] = fminunc(fun1,x0);
d = xmin1(5);

global d;global t; global U;
init=zeros(4,48);
x0 = xmin1(1:4,1);

fun2 = @ss_init2;
[xmin2,fmin2,flag2] = fminunc(fun2,x0);

for t=1:48
init(:,49-t)=xmin2;
end

xi0(:,1)=reshape(init,192,1);

global T;global Y;global xi0; global xi1; global P0; global P1; global F;
global d11;global d12;global d13;global Q;global R;global D;

% function
fun = @ekf_chv_baseline % used for the baseline model (estimated with baseline/alternative/unadjusted/post94 data)
option = optimset('MaxIter',1000000,'MaxFunEvals',1000000,'TolX',1e-4)

% starting values %
% alternative dataset :Hornstein (2012)'s
 x0 = [-0.6279,-0.614,-0.6119,-0.791,2.9993,-0.1531,1.7208,-0.2206,-0.2625,-0.27,-0.3308,-0.3778];
[xmin_u1,fval_u,exitFlag_u,output_u] = fminunc(fun,x0,option);

% smoothed inference %%
[smW,smPU,sxi] = smoothin_base(F,P0,P1,xi0,xi1,T);

%% Standard error calculation
% hessian and jacobian for the baseline model
program1 = @ekf_chv_hes; % hessian
program2 = @ekf_chv_jac; % jacobian

sqs1 = sqrt(Q(1,1));sqs2 = sqrt(Q(2,2));sqs3 = sqrt(Q(3,3));sqs4 = sqrt(Q(4,4));
sqr1 = sqrt(R(1,1));sqr2 = sqrt(R(2,2));sqr3 = sqrt(R(3,3));sqr4 = sqrt(R(4,4));sqr5 = sqrt(R(5,5));
c1=d11;c2=d12;c3=d13;
x0 = [sqs1,sqs2,sqs3,sqs4,c1,c2,c3,sqr1,sqr2,sqr3,sqr4,sqr5]'; 

[hess,err]=hessian(program1,x0);
[white_SE]=white(program2,x0,hess);   

%% Variance Decomposition

[var_share] = vardcmp(Q,sxi,D); % baseline/ Hornstein/ unadjusted
var_q = var_share;
var_io = [var_share(:,1)+var_share(:,2),var_share(:,3)+var_share(:,4)];
var_HL = [var_share(:,1)+var_share(:,3),var_share(:,2)+var_share(:,4)];

program = @ekf_chv_base_sim

%% Standard errors for variance decomposition
nit=1000; % number of iteration

mle = x0;
sigma = inv(-hess);
sigma2 = (sigma + sigma.')/2; 
rng default  % For reproducibility
rv = mvnrnd(mle,sigma2,nit);

global Y;global T;global xi0;

for i=1:nit
    i
x0 = rv(i,:);
tic
[rq3,rq12,rq24]=feval(program,x0);%function evaluation
toc
 Q_3m(i,:) = rq3;
 Q_12m(i,:) = rq12;
 Q_24m(i,:) = rq24;
end

table4_se(1:4,1) = std(Q_3m,0,1)';
table4_se(5:8,1) = std(Q_12m,0,1)';
table4_se(9:12,1) = std(Q_24m,0,1)';


save result_hs_3prm_mat