% minPosterior.m

addpath functions

%% starting values for parameters

x0_ = [0.02*10      % sig_g 
       0.0015*10    % sig_tauc
       0.005*10     % sig_tauw
       0.008*10     % sig_tauK
       0.1*10       % sig_tau 
       0.002*10     % sig_z
       0.003*10     % sig_rstar
       0.02*10      % sig_zeta
       0.02*10      % sig_L
       0.08*10      % sig_I
       0.01*10      % sig_b
       3            % Upsilon
       0.15         % h
       0.1          % psi_b
       0.05*10      % gam_gb
       0.05*10      % gam_cb
       .1*10        % gam_wb
       0.01*10      % gam_Kb
       0.03*10      % gam_taub
       0.3          % rho_g
       0.5          % rho2_g
       0.3          % rho_c
       0.5          % rho2_c
       0.3          % rho_w
       0.5          % rho2_w
       0.3          % rho_K
       0.5          % rho2_K
       0.3          % rho_tau
       0.5          % rho2_tau
       0.3          % rho_z
       0.7          % rho_r_star
       0.9          % rho_zeta
       0.95         % rho_L
       0.25         % rho_I
       0.36         % rho_b
       -.0035*100   % Delta z
       2.1          % mu_sig
       -0.2         % Delta_tauw
       0.2          % Delta_tauK
       ];


dim_x = size(x0_,1);

param_ = param_func(x0_);

S_.Options.EST.kron = 0;   % Set to 1 for Sigma = (eye((n_)^2)-(kron(Pt(:,:,1),Pt(:,:,1))))\DOmega(:);

%% Use FMINCON

warning('off','all') ; 

aus_data;  
S_.data = Zobs; 

post_ = @(x_,T_) loglike(x_,T_,S_) + prior(x_);

post_min = @(x) post_(x,[88, 83]);

% Add non-negativity constraints.
A =[];
b = [];
Aeq = [];
beq = [];
lb = zeros(dim_x,1);
lb(param_.order.Delta_z) = -1;
lb(param_.order.Delta_tauw) = -3; 

ub = inf*ones(dim_x,1);
ub(param_.order.rho_g:param_.order.rho_I) = 0.99999999;


% Options
options = optimset;
% Modify options setting
options = optimset(options,'Display', 'off');
options = optimset(options,'TolFun', 1e-4);
options = optimset(options,'TolX',1e-6); 
options = optimset(options,'PlotFcns', {  @optimplotfunccount @optimplotfval @optimplotstepsize @optimplotx});
options = optimset(options,'Algorithm', 'interior-point');
options = optimset(options,'Diagnostics', 'off');
options = optimset(options,'MaxIter',100);


[x_mode, post_val, exitflag,output,lambda,grad,hess_] = fmincon(post_min,x0_,A,b,Aeq,beq,lb,ub,[],options);

save('hessmat','hess_','x_mode','post_val');