clear variables
global Kss Lss Yss Css Rdss Riss Rss debtss wss r_Kss ...
       Pi_ss Pi_starss Ass SGss Sss vss qss Iss x1ss x2ss mcss

global R_star betta siggma phi theta ksi psi_pi psi_y sg rhoA rhoG ...
    chi_d rhoS dbar Rbar chi_R delta alppha chi_I mu rho

dbar = 0.1438373; 
Rbar = 0.92;
Pi_ss = 1.01;

%%%%%%%%%%%% Don't mind this first part, is just data preparation %%%%%%%%%


%%%%%% Arranging Data - input is the soe_nk_data and trend coefficients
load data
varobs = 5;
nobs = 83;
trends_coefs = 3; %quadratic trends, including constant

trends = zeros(varobs, trends_coefs);
trends(1,:) = [2.071168e-01  4.650618e-03  4.783001e-05 ]; %R_Y trend
trends(2,:) = [ 2.284219e-01 -6.848424e-03  7.624129e-05 ]; %d_Y trend
trends(3,:) = [ 7.824197431   0.019958475  -0.000120018 ]; %logY trend
trends(4,:) = [5.9941129251  0.0300870340 -0.0002582194]; %logInv trend
trends(5,:) = [1.020079e+00 -1.153809e-04 -7.474502e-08]; %Pi trend

% Adjut R_Y, d_Y and Pi so they fluctuate around their steady state
% Adjust logY and logInv so that they fluctuate around 0
time_index = 1:nobs;

%Save it to use later
RY_obs_trended = RY_obs;
dY_obs_trended = dY_obs;
logY_obs_trended = logY_obs;
Pi_obs_trended = Pi_obs;

for i = 1:nobs
    RY_obs(i) = RY_obs(i) - (trends(1, 1) + time_index(i)*trends(1, 2) + trends(1, 3)*time_index(i)^2); %+ Rbar;
    dY_obs(i) = dY_obs(i) - (trends(2, 1) + time_index(i)*trends(2, 2) + trends(2, 3)*time_index(i)^2); %+ dbar;
    logY_obs(i) = logY_obs(i) - (trends(3, 1) + time_index(i)*trends(3, 2) + trends(3, 3)*time_index(i)^2);
    logInv_obs(i) = logInv_obs(i) - (trends(4, 1) + time_index(i)*trends(4, 2) + trends(4, 3)*time_index(i)^2);
    Pi_obs(i) = Pi_obs(i) - (trends(5, 1) + time_index(i)*trends(5, 2) + trends(5, 3)*time_index(i)^2); %+ Pi_ss;
end

% Save datafile
save soe_nk_data RY_obs dY_obs logY_obs logInv_obs Pi_obs

%%%% Declaring parameters. 
%Initial values are the mean of posteriors of previous exercise

R_star = 1.003742;
betta = 1/R_star;
siggma = 2;
phi = 1.3077;
theta = 0.6413;
ksi = 6.1394;
psi_pi = 2.5413;
psi_y = 0.2167;
sg = 0.1948;
rhoA = 0.9319 ; sigmaA = 0.0096;
rhoG = 0.8340; sigmaG = 0.2437;
%Pi_ss = 1.01;
chi_d = 0.5595;
rhoS = 0.8361; sigmaS = 0.0046;
%dbar = 0.1438373; 
%Rbar = 0.8276067; 
chi_R =  0.0274;
delta = 0.015;
alppha = 0.41;
chi_I = 2.6019;
rho = 0.5139;

save parametersfile  R_star betta siggma phi theta ksi psi_pi psi_y sg rhoA sigmaA ...
    rhoG sigmaG Pi_ss chi_d rhoS sigmaS dbar Rbar chi_R chi_I delta alppha mu rho
    

%%%% Setting Steady State of the model
Riss = Pi_ss/betta;
Pi_starss = ((Pi_ss^(1-ksi) - theta)/(1-theta))^(1/(1-ksi));
vss = (1-theta)/(1-theta*Pi_ss^ksi)*(Pi_ss/Pi_starss)^ksi;
mcss = (1-betta*theta*Pi_ss^ksi)/(1-theta*betta*Pi_ss^(ksi-1))*Pi_starss/Pi_ss*(ksi-1)/ksi;
qss = 1;
r_Kss = 1/betta + delta -1;
Rdss = R_star;
wss = (1-alppha)*mcss^(1/(1-alppha))*(alppha/r_Kss)^(alppha/(1-alppha));
R_Y = Rbar; d_Y = dbar;
K_L = alppha/(1-alppha)*wss/r_Kss;
Aux = delta*K_L + ( (1-R_star)*R_Y + (Rdss-1)*d_Y - (1-sg) )*K_L^alppha*1/vss;
Lss = (wss/((-Aux)^siggma))^(1/(phi+siggma));
Kss = K_L*Lss;
Iss = delta*Kss;
Yss = Kss^alppha*Lss^(1-alppha)*1/vss;
Css = -Aux*Lss;
Rss = R_Y*Yss;
debtss = d_Y*Yss;
x1ss = Css^(-siggma)*mcss*Yss/(1 - betta*theta*Pi_ss^ksi);
x2ss = ksi/(ksi-1)*Pi_ss/Pi_starss*x1ss;
Ass = 1; Sss= 1; SGss = sg;

%check_ss = steady_check(); %steady state seems correct. 

path(path,'C:\dynare\4.6.3\matlab')
dynare SOE_NK_estimation3;



% Thanks for your answer ! You are right, higher order solutions are implying explosives paths when the shock is "big" enough, but solutions are fine when I use the pruning option in the stoch_simul command. However, I still get the same problem after setting :
% 
% > options_.particle.pruning=1;
% 
% After ***varobs*** and before ***estimation*** command. This is strange since solutions appears to be fine after pruning process. To see just execute the .m file ***main_NK_estimation3.m***. Since first-order works fine I am tempting to let go the estimation for higher orders solutions. Rather I might use higher-order solutions with parameters calibrated from this first-order exercise. 
% 
% Seizing the opportunity here, is there any way to generate the output of forecasts ***oo_.MeanForecast*** 
% 


