clear all 

%add Dynare folder 
%addpath('C:\Users\JP\Dropbox\Gemeinsame Dateien Benni\dynare\4.4.3\matlab')

%run Dynare once to set up files and structures
dynare Model_Uncertainty_Moment1 noclearall
%keep only important variables in workspace
keep('M_','options_','oo_')

%get variable indices for moment matching
vars={'y_obs';'c_obs';'pi_obs';'ffr_obs';'h_obs'};
varsindex=[];
for ii = 1:size(vars,1)
  varname = deblank(vars(ii,:));
  for jj=1:M_.orig_endo_nbr
     if strcmp(varname,deblank(M_.endo_names(jj,:)))
         varsindex=[varsindex; jj];
     end
  end
end

%load data
load('data_first_diff_BP')

% compute moments for matching
empirical_data_cycle = [y_obs c_obs pi_obs ffr_obs h_obs];
empiricalmoments(:,1)=std(empirical_data_cycle);
empiricalmoments(:,2)=std(empirical_data_cycle)./std(empirical_data_cycle(:,1));
empiricalmoments(:,3)=corr(empirical_data_cycle,empirical_data_cycle(:,1));
empiricalmoments(:,4)=diag(corr(empirical_data_cycle(2:end,:),empirical_data_cycle(1:end-1,:)));
empiricalmoments_table=empiricalmoments;

maxq = floor(4*(size(empirical_data_cycle,1)/100)^(2/9));
optwt = 1; % 2 if optimal, 1 if diagonal, 0 if identity
[empiricalmoments,wtmat] = Wtmat(empirical_data_cycle,optwt,maxq); % data moments and weighting matrix

%define length of simulation and other options
options_.long=10*length(empirical_data_cycle);
options_.drop=1500;
options_.pruning=1;

%% generate random shocks
ii=1;
smmstream = RandStream('mt19937ar','Seed',ii);
shocks=randn(smmstream,options_.long+options_.drop,M_.exo_nbr);

%set starting values for parameters
h=0.81;                                                                                    
rho_er=0.81; 
rho_barpi=0.995;
alpha=0.6;                                                                                          
rho=0.5;                                                                                                                                            
phi_pi=1.7768;                                               
phi_y=0.31949;                                                                                                                                          

%transform to unbounded domain
par=partransform_unbounded([h;rho_er;rho_barpi;alpha;rho;phi_pi;phi_y]');
                           

%set options for csminwel
H0 = 1e-4*eye(length(par));
crit = 1e-8;
nit = 1000;
verbose = 2;

%run csminwel
[fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
        csminwel(@smm_objective_function,par,H0,[],crit,nit,empiricalmoments,varsindex,shocks,M_,options_,oo_,wtmat);

save results
