

global M_ oo_ options_

% Fetch matrix of shocks
load 'shockmatrix_for_results'; %Use this matrix to reproduce estimates

shocks=shockmatrix;

Simulnumber=1;
nreplic = Simulnumber;

Simulperiods=length(shocks);
Discardperiods=0;

M=M_; oo=oo_; options=options_;

TShadow=20;
TIRF=20;
load muIss; load muEss;
steadymuvalues=[ muIss muEss ];
var_list_=[];

Model7Data=StartDynSimM7(TShadow,TIRF,'muihat','muehat',steadymuvalues,var_list_ );

ZeroVec = Model7Data.ZeroVec;
OneVec = Model7Data.OneVec;
AlphaStart = Model7Data.AlphaStart;
Options = Model7Data.Options;
Var_cons1 = Model7Data.Var_cons1;
Var_cons2 = Model7Data.Var_cons2;
Accuracy = Model7Data.Accuracy;
TshadowI = Model7Data.TshadowI;
Mmat = Model7Data.Mmat( 1:TshadowI*2, : ); % only use the top block (Should be redundant if TshadowI=TIRF) 
vars = Model7Data.Vars;
sDM = Model7Data.sDM;


    
%Stats=zeros(nreplic,65);
SimOutputUnfiltered = zeros(size( oo.dr.ys, 1 ),Simulnumber,Simulperiods-Discardperiods);


jj=nreplic;


    
TSimul=Simulperiods;


    
% Generate vectors of innovations
    
nper=TSimul+1;



 	tmp_exo_simul = zeros( TSimul, M.exo_nbr );
 	oo.exo_simul = tmp_exo_simul;
    oo.exo_simul(:,1:3)=shocks';

    
	shadow_indexI = strmatch( 'shadow_sumImp_0', M.endo_names, 'exact' ); 
	epsilon_shadowI_index = strmatch( 'epsilon_shadowImp_0', M.exo_names, 'exact' );
        shadow_indexE = strmatch( 'shadow_sumEnt_0', M.endo_names, 'exact' ); 
	epsilon_shadowE_index = strmatch( 'epsilon_shadowEnt_0', M.exo_names, 'exact' );

	ex_ = oo.exo_simul;
	ex_( :, epsilon_shadowI_index:end ) = 0;
	dr = oo.dr;
	iter = size( ex_, 1 );
	k2 = dr.kstate( dr.kstate( :, 2 ) <= M.maximum_lag+1, [ 1 2 ] );
	k2 = k2( :, 1 )+( M.maximum_lag+1-k2( :, 2 ) )*M.endo_nbr;
	order_var = dr.order_var;

	Var_cons1_index = strmatch( Var_cons1, M.endo_names, 'exact' );
        Var_cons1_steady = muIss;
        Var_cons2_index = strmatch( Var_cons2, M.endo_names, 'exact' );
        Var_cons2_steady = muEss;
	SIMU_CONS = zeros( size( dr.ys, 1 ), iter+M.maximum_lag );
	SIMU = SIMU_CONS;
	
    usize=size( dr.ghu );
    save('usize','usize');
	oghu = zeros( usize );
    save('order_var','order_var');
	oghu( order_var, : ) = dr.ghu;
	ZIMat = zeros( 2*TshadowI, M.exo_nbr );
        ZIMat( 1:TshadowI, ( end - TshadowI*2 + 1 ):end-TshadowI ) = eye( TshadowI );
	ZIMat( TshadowI+1:end, ( end - TshadowI + 1 ):end ) = eye( TshadowI );
	if any( abs( oghu( shadow_indexI:end, : ) - ZIMat ) > Accuracy * sqrt( Accuracy ) )
		error( 'Model3:IncorrectShadowSpecification', 'You have specified shadow shocks incorrectly.' );
	end
	oghu( shadow_indexI:end, : ) = ZIMat;

	% Starting the simulation

	disp( 'Starting simulation. Please be patient.' );

	LongestBindPeriod = 0;
	CurrentBindPeriod = 0;
		for ii = 2:iter+M.maximum_lag
		yhat_cons = SIMU_CONS( order_var( k2 ), ii-1 );
		yhat = SIMU( order_var( k2 ), ii-1 );
		
		lag_term = zeros( length( dr.ys ), 1 );
		lag_term( order_var ) = dr.ghx * yhat;
		
		lag_term_cons = zeros( length( dr.ys ), 1 );
		lag_term_cons( order_var ) = dr.ghx * yhat_cons;

		epsilon = ex_( ii-1, : )';
		epstmp = oghu * epsilon;
		
		SIMU( :, ii ) = lag_term + epstmp;
		
		epsilon( epsilon_shadowI_index:end ) = -lag_term_cons( shadow_indexI:end );
		temp_current_SIMU_CONS = lag_term_cons + oghu * epsilon;
        
        
            
            
        Vi=zeros(TshadowI,1);
        Ve=zeros(TshadowI,1);
            
        Vi(1) = Var_cons1_steady + (Var_cons1_steady.*temp_current_SIMU_CONS( Var_cons1_index ));
        Ve(1) = Var_cons2_steady + (Var_cons2_steady.*temp_current_SIMU_CONS( Var_cons2_index ));
        
        for j=2:TshadowI
            temp_current_SIMU_CONS( order_var ) = dr.ghx * temp_current_SIMU_CONS( order_var( k2 ) );
	    Vi( j ) = Var_cons1_steady + (Var_cons1_steady.*temp_current_SIMU_CONS(Var_cons1_index ));
            Ve( j ) = Var_cons2_steady + (Var_cons2_steady.*temp_current_SIMU_CONS( Var_cons2_index ));
        end
        
        V = [Vi; Ve];
        

		alpha = solve_quadprog_both( Mmat, V, AlphaStart, sDM, TshadowI, Accuracy, Options, ZeroVec, OneVec );

		epsilon( epsilon_shadowI_index:end ) = epsilon( epsilon_shadowI_index:end ) + alpha;

		SIMU_CONS( :, ii ) = lag_term_cons + oghu * epsilon;

		
		if alpha( 1 ) > Accuracy
			CurrentBindPeriod = CurrentBindPeriod + 1;
		else
			if CurrentBindPeriod > LongestBindPeriod
				LongestBindPeriod = CurrentBindPeriod;
			end
			CurrentBindPeriod = 0;
        end
        save('ii','ii');
	end

	if LongestBindPeriod > TshadowI
		warning( 'Model7:Inaccuracy', 'Simulations may be inaccurate, as the constraint may be non-binding for more than TshadowI periods along the simulated path.' );
	end

	% Save results
    SIMU = bsxfun( @plus, SIMU, dr.ys );
    SIMU_CONS = bsxfun( @plus, SIMU_CONS, dr.ys );
    % The ONBC simulation is in SIMU_CONS
    
    % Transform shadow prices to level values, so as to visually inspect
    % that ONBC is working and shadow prices are bounded at zero
    SIMU_CONS(8,:)=(muIss.*SIMU_CONS(8,:))+muIss;
    SIMU_CONS(9,:)=(muEss.*SIMU_CONS(9,:))+muEss;
    % Do the same transformation with the unbounded simulation to compare
    SIMU(8,:)=(muIss.*SIMU(8,:))+muIss;
    SIMU(9,:)=(muEss.*SIMU(9,:))+muEss;
    
    SIMU=SIMU(:,2:end);
    SIMU_CONS=SIMU_CONS(:,2:end);

	for varIndex = 1:size( vars, 1 )
		var = strtrim( vars( varIndex, : ) );
		SIMU_.( var ) = SIMU( strmatch( var , M.endo_names, 'exact' ), 2:end );
		SIMU_CONS_.( var ) = SIMU_CONS( strmatch( var , M.endo_names, 'exact' ), 2:end );
    end
    

SIMU_CONS=SIMU_CONS(:,(Discardperiods+1):end);
SIMU=SIMU(:,(Discardperiods+1):end);

perids=length(SIMU_CONS(1,:));

indic_i=zeros(1,perids);
indic_e=zeros(1,perids);
indic_both=zeros(1,perids);

for pp=1:perids
    if SIMU_CONS(8,pp)<0.00000000001
        indic_i(1,pp)=1;
    end
    if SIMU_CONS(9,pp)<0.00000000001
        indic_e(1,pp)=1;
    end
    if (SIMU_CONS(8,pp)<0.00000000001) && (SIMU_CONS(9,pp)<0.00000000001)
        indic_both(1,pp)=1;
    end
end


SimOutputUnfiltered(1:length(SIMU_CONS(:,1)),jj,:) = SIMU_CONS;


save('SIMU_CONS','SIMU_CONS');

% Produce moments to be matched

vv= [ SIMU_CONS(24,:)' SIMU_CONS(28,:)' SIMU_CONS(23,:)' SIMU_CONS(18,:)' SIMU_CONS(30,:)'  ];
vv = [vv(5:end,1:4)-vv(5-4:end-4,1:4) vv(5:end,5)];
          
save('vv','vv');

SimLevels = [ SIMU_CONS(24,:); SIMU_CONS(28,:); SIMU_CONS(23,:); ];
save('SimLevels','SimLevels');
