%% run the factor model with the baseline model %%
clear all;
load Data;

% use the appropriate data: sample period, (Feb 1972 - June 2017)
U = U_baseline; % the baseline 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;
global l1; global l2; global l3; global l4;
global s1; global s2; global s3; global s4;

% function
fun = @ekf_chv_fc % used for the factor model
option = optimset('MaxIter',1000000,'MaxFunEvals',1000000,'TolX',1e-4)

% starting values : factor model
x0 = [-0.0032,-0.0028,0.0085,0.0008,-0.6575,-0.6453,-6.0812,-0.7975,4.7401,-0.0496,2.0382,-0.217,-0.2795,-0.2656,-0.3323,-0.4094];
[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 factor model
 program1 = @ekf_chv_fc_hes; % hessian
 program2 = @ekf_chv_fc_jac; % jacobian
 
x=xmin_u1;
l1=-1+2*exp(10*x(1))/(1+exp(10*x(1)));
l2=-1+2*exp(10*x(2))/(1+exp(10*x(2)));
l3=-1+2*exp(10*x(3))/(1+exp(10*x(3))); 
l4=-1+2*exp(10*x(4))/(1+exp(10*x(4)));
L = [l1;l2;l3;l4];LL = L*L';
s1=exp(10*x(5))/(1+exp(10*x(5)))+0.000001;
s2=exp(10*x(6))/(1+exp(10*x(6)))+0.000001;
s3=exp(10*x(7))/(1+exp(10*x(7)))+0.000001; 
s4=exp(10*x(8))/(1+exp(10*x(8)))+0.000001;
ss = [s1,0,0,0;...
       0,s2,0,0;...
       0,0,s3,0;...
       0,0,0,s4];
ss1 = sqrt(s1);ss2 = sqrt(s2);ss3 = sqrt(s3); ss4 = sqrt(s4);
c1=d11;c2=d12;c3=d13;
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));
x0 = [l1,l2,l3,l4,ss1,ss2,ss3,ss4,c1,c2,c3,sqr1,sqr2,sqr3,sqr4,sqr5];

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

%% Variance Decomposition

[var_share] = vardcmp_fc(LL,ss,sxi,D);
fvar_q = var_share;
fvar_io = [var_share(:,1),var_share(:,2)+var_share(:,3),var_share(:,4)+var_share(:,5)];
fvar_HL = [var_share(:,1),var_share(:,2)+var_share(:,4),var_share(:,3)+var_share(:,4)];

%% Standard errors for the variance decomposition

program = @ekf_chv_fc_sim

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:5,1) = std(Q_3m,0,1)';
table4_se(6:10,1) = std(Q_12m,0,1)';
table4_se(11:15,1) = std(Q_24m,0,1)';

 save result_base_3prm_fc.mat