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
   
%%%% 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.2\matlab')
dynare SOE_NK_estimation2;

%%%% Ignore this below

% %%%%%%% Forecasts
% %%% Reserves
% nobs = 83; %number of observations
% nforecast = 40; %number of forecast periods ahead
% RY_forecast = zeros(nforecast, 3); 
% RY_forecast(:, 1) = oo_.MeanForecast.HPDinf.RY_obs;
% RY_forecast(:, 3) = oo_.MeanForecast.HPDsup.RY_obs;
% RY_forecast(:, 2) = oo_.MeanForecast.Mean.RY_obs;
% 
% trendR = zeros(1, 3); %recover quadratic trend coefficients
% trendR(1, 1) = 2.126026e-01; %const
% trendR(1, 2) = 4.101705e-03; %t
% trend(1, 3) = 5.775059e-05; %t^2
% t_basis = nobs + 1; %t to include back the trend
% 
% RY_forecast = RY_forecast - Rbar; %remove steady state constant
% 
% for i  = 1:nforecast
%     RY_forecast(i, :) = RY_forecast(i, :) + trendR(1,1) + trendR(1,2)*t_basis + trendR(1,3)*t_basis^2;
%     t_basis = t_basis + 1;
% end
% 
% % Ploting results
% alpha = 0.5; %gray scale, varies between [0,1]
% periods = 1:40;
% plot(periods, RY_forecast(:, 1), 'color', [0,0,0] + alpha, 'LineWidth', 2)
% hold on
% plot(periods, RY_forecast(:, 2), 'b', 'LineWidth', 2)
% hold on
% plot(periods, RY_forecast(:, 3), 'color', [0,0,0] + alpha, 'LineWidth', 2)
% 

