function objective_fun=smm_objective_function(par,momentstomatch,varsindex,shocks,M_,options_,oo_,wtmat)
% function outvalue=smm_objective_function(par,momentstomatch,varsindex,shocks,M_,options_,oo_,wtmat)
% Inputs_
% - par             [vector]    vector of parameters to estimate
% - momentstomatch  [vector]    vector of target moments
% - varsindex       [vector]    vector of indices for variables in Dynare
%                               decision rules/simulations
% - shocks          [matrix]    matrix of random numbers for simulations
% - M_              [structure] Dynare model structure
% - options_        [structure] Dynare options structure
% - oo_             [structure] Dynare results structure
% - wtmat           [matrix]    weighting matrix in objective
%
% Outputs:
% - objective_fun   [scalar]    value of the objective function
%
% written by Johannes Pfeifer and Benjamin Born

%retransform parameters to bounded support
h=1/(1+exp(par(1)));
rho_er=1/(1+exp(par(2)));
rho_barpi=1/(1+exp(par(3)));
alpha    = 1/(1+exp(par(4)));
rho       = 1/(1+exp(par(5)));
phi_pi     = 1.5+exp(par(6));
phi_y      = exp(par(7));




%save them in vector for prior computation later on
theta=[h;rho_er;rho_barpi;alpha;rho;phi_pi;phi_y];

%set parameters for simulation in Dynare
set_param_value('h',h)

set_param_value('rho_er',rho_er)
set_param_value('rho_barpi',rho_barpi)
set_param_value('alpha',alpha)
set_param_value('rho',rho)

set_param_value('phi_pi',phi_pi)
set_param_value('phi_y',phi_y)


tic

%get decision rules
if isempty(options_.qz_criterium)
    options_.qz_criterium = 1+1e-6; %make sure option is set
end
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
oo_.dr=dr;

if info==0 %did solve
    %% simulate series
    temps = repmat(oo_.dr.ys,1,M_.maximum_lag); %inital state vector
    y1 = simult_(temps,oo_.dr,shocks,options_.order); %simulate

    if any(any(isnan(y1))) || any(any(isinf(y1))) %simulation exploded, should never happen
        objective_fun= 10^10
        return
    end

    simulated_data_cycle = y1(varsindex,end-options_.long:end)'; %get the simulated values

%     %compute the model moments
%     modelmoments(:,1)=std(simulated_data_cycle);
%     modelmoments(:,2)=std(simulated_data_cycle)./std(simulated_data_cycle(:,1));
%     modelmoments(:,3)=corr(simulated_data_cycle,simulated_data_cycle(:,1));
%     modelmoments(:,4)=diag(corr(simulated_data_cycle(2:end,:),simulated_data_cycle(1:end-1,:)));

    % Variance and covariances
    nvar=length(varsindex);
    vars = simulated_data_cycle(3:options_.long,end);
    vars = [vars simulated_data_cycle(3:options_.long,:).*simulated_data_cycle(3:options_.long,:)];

    for jj = 2:nvar %covariance with output
         vars = [vars simulated_data_cycle(3:options_.long,1).*simulated_data_cycle(3:options_.long,jj)];  
    end

    %first and second order autocovariance
    autoc1 = simulated_data_cycle(3:options_.long,:).*simulated_data_cycle(2:options_.long-1,:);
    autoc2 = simulated_data_cycle(3:options_.long,:).*simulated_data_cycle(1:options_.long-2,:);

    moments_model= [mean(vars)';mean(autoc1)';mean(autoc2)'];
    
    %compute objective function
    moments_difference=moments_model-momentstomatch;
    objective_fun=moments_difference'/wtmat*moments_difference+compute_prior_weights(theta)

    if isnan(objective_fun) || abs(objective_fun)==Inf
        objective_fun=10^6+compute_prior_weights(theta);
    end
    toc
else %did not solve  
    objective_fun= 10^6+compute_prior_weights(theta);
end
