
format short

%restoredefaultpath
addpath 'D:\_1_PhD_Research\OBC_Iacoviello\occbin_20140630\Occbin_update-4.5\toolkit_files'
%setpathdynare4

clear
set(0,'DefaultLineLineWidth',2)
close all



load PARAM_EXTRA_CALIBRATED
load mle_estimates_fminsearch params_labels params* filtered_errs sample_length err_list
for i=1:numel(params_labels)
    evalc([ cell2mat(params_labels(i)) '= params1(' num2str(i) ')']) ;
end

XX = [ NaN ]
save  PARAM_EXTRA XX % these two rows no use



    
    
    
    STD_J      =0.0737;
    STD_K      =0.0360;
    STD_P      =0.0030;
    STD_R      =0.0013;
    STD_W      =0.0100;
    STD_Z      =0.0163;
    
    
    
    modnam_00 = 'baby00'; % base model (constraint 1 and 2 below don't bind)
    modnam_10 = 'baby10'; % first constraint is true
    modnam_01 = 'baby01'; % second constraint is true
    modnam_11 = 'baby11'; % both constraints bind
    
    % compile constraint
    constraint1 = 'lm<-lm_ss';
    constraint_relax1 = 'maxlev > 0';
    constraint2 = 'rnot < -log(r_ss)';
    constraint_relax2 = 'r > -log(r_ss)';
    
    
    curb_retrench = 0;           % if 1 slow down relaxation of constraints
    maxiter = 10;                % number of iterations allowed to look for a solution
    
    %irfshock =char('eps_j','eps_k','eps_p','eps_r','eps_w','eps_z','eps_j4');
    
    
    
    %sequence1 = [ STD_J STD_K STD_P STD_R STD_W STD_Z ];
    %nperiods = 40;
    
    
    irfshock =char('eps_j4','eps_j');
 
%initial_condition_states = repmat(oo_.dr.ys,1,M_.maximum_lag);
%shock_matrix = zeros(options_.irf,M_.exo_nbr); %create shock matrix with number of time periods in columns
%shock_matrix(2,strmatch('eps_j4',M_.exo_names,'exact')) = 0.1; %set news shock to 1 (use any shock size you want)
%shock_matrix(2+4,strmatch('eps_j',M_.exo_names,'exact')) = -0.1; %8 periods later use counteracting shock of -1
shock_matrix = zeros(40,2);
shock_matrix(2,1) = 0.1;
shock_matrix(2+4,2) = -0.1;
shockssequence= shock_matrix;
   nperiods = numel(shockssequence);    
     
    
    [zdatal zdatap zdatass oobase_ Mbase_] = solve_two_constraints(...
        modnam_00,modnam_10,modnam_01,modnam_11,...
        constraint1, constraint2,...
        constraint_relax1, constraint_relax2,...
        shockssequence,irfshock,nperiods+50,curb_retrench,maxiter);
    
    
    
    
    for i=1:Mbase_.endo_nbr
        eval([deblank(Mbase_.endo_names(i,:)),'_l=zdatal(1:nperiods,i);']);
        eval([deblank(Mbase_.endo_names(i,:)),'_p=zdatap(1:nperiods,i);']);
        eval([deblank(Mbase_.endo_names(i,:)),'_ss=zdatass(i);']);
    end
    
    for i = 1:size(Mbase_.param_names,1)
        eval([Mbase_.param_names(i,:),'= Mbase_.params(i);']);
    end
 for i = 1:size(Mbase_.param_names,1)
        eval([Mbase_.param_names(i,:),'= Mbase_.params(i);']);
    end
    
    ctot_ss = c_ss + c1_ss;
    ctot_l = c_l + c1_l;
    ctot_p = c_p + c1_p;
    lev_p = (exp(b_p/b_ss)*b_ss)./(exp(q_p/q_ss)*q_ss)./(exp(h1_p/h1_ss)*h1_ss);
    levagg_p = (exp(b_p/b_ss)*b_ss)./(exp(q_p/q_ss)*q_ss);
    lev_l = (exp(b_l/b_ss)*b_ss)./(exp(q_l/q_ss)*q_ss)./(exp(h1_l/h1_ss)*h1_ss);
    
    
    titlelist = char(...
        'Interest Rate',...
        'GDP',...
        'House Prices',...
        'Total Consumption',...
        'Consumption Saver',...
        'Consumption Borrowers',...
        'Investment',...
        'Multiplier (level)',...
        'Hours');
    ylabels = char('Percentage Point, AR','Percent','Percent',...
        'Percent','Percent','Percent',...
        'Percent','Level','Percent');
    
    
subplot(2,2,1)
plot(y_p*100/y_ss,'b','LineWidth',2);  hold on
title('output')
ylabel('%from steady state');  % if you add %from steady state, you should*100
hold on;

subplot(2,2,2)
plot( q_p*100/q_ss,'b','LineWidth',2);  hold on
title('House Price')
ylabel('%from steady state');
hold on;

    
