%% A1 : estimation %% 

clear all;

addpath('../data')

load ne_data.mat

Y= U(48:end,:)/1000000;

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

%% initial vector computation : already calculated in the dataset.
%{
global T;global Y;global nt;
nt = 48;
x0 = [0.2;1;-2;0.1;0.01];
lb = [-5;-5;-5;-5;-0.5];
ub =  [5;10;5;5;0.5];
[xmin1,fmin1,flag1] = fminunc(@ss_ph48_all,x0);
% Y = Y(1:end-2,:);
% [xmin1,fmin1,flag1] = fmincon(@ss_ph48_all,x0,[],[],[],[],lb,ub);
% [xmin1,fmin1,flag1] = fminunc(@ss_ph48_all,x0);

d = xmin1(5);

global d;global t; global U;
U=U/1000000;
init=zeros(4,48);
x0 = [0.1;1;-2;0.1];
lb = [-5;-5;-5;-5];
ub =  [5;10;5;5];
[xmin2,fmin2,flag2] = fmincon(@ss_ph48_tv,x0,[],[],[],[],lb,ub);
for t=1:48
init(:,49-t)=xmin2;
end

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

%}

%% estimation of nonlinear state space model

global T;global Y;global xi0; global xi1; global P0; global P1; global F;
global d11;global d12;global d13;global d14;global Q;global R;global BD;global RD;
global r11;global r12;global r13;global r14;global ind65;

load indfile
ind65 = ind2(size(U,1)-T+1:size(U,1),1);
fun = @ekf_PH3d_rec
option = optimset('MaxIter',1000000,'MaxFunEvals',1000000,'TolX',1e-6)

x0 = [-1.1299;-0.6576;0.0403;-0.6664; ...
 -0.6052;0.2356;0.4354;-0.4244;-0.4552;
 -0.5618;-0.6068;-0.7189;-0.6145;1.7517;-0.5553];

[xmin_u2,fval_u2,exitFlag_u2,output_u1] = fminunc(fun,x0,option);

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

% variance decomposition
[var_share] = vardcmp_rec(Q,sxi,BD,RD);
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)];

save('../result_files/ne_3prm_rec_e24_final.mat')

%% type L probability and path of genuine duration dependence.

%{
for i=1:48
for t=1:409 % end point can change 
    if ind65(t,1)==0
    D=BD;
    else
    D=RD;
    end
smp(t,i)=exp(-exp(sxi(3,t)+BD(i,1)));
end
end
%}