
tic;
nreplic = Simulnumber;

global M_
global oo_

SimOutputLevels = zeros(73,Simulnumber,Simulperiods-Discardperiods);

    TSimul=Simulperiods;

	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;
	M = Model7Data.M( 1:TshadowI*2, : ); % only use the top block (Should be redundant if TshadowI=TIRF)
	vars = Model7Data.Vars;
	sDM = Model7Data.sDM;
    dr = oo_.dr;
    
    endo_names  = M_.endo_names; % To avoid a global variable in parfor loop
    exo_names   = M_.exo_names ;
    endo_nbr    = M_.endo_nbr  ;
    exo_nbr     = M_.exo_nbr   ;
    maximum_lag = M_.maximum_lag ;


parfor jj=1:nreplic 
    
% Generate vectors of innovations
    
nper=TSimul+1;


ea=randn(1,nper)*stdev_ea;
ee=randn(1,nper)*stdev_eps;
es=randn(1,nper)*stdev_es;

eps = [ es ; ea ; ee ];
eps = eps(:,2:end);


 	tmp_exo_simul = zeros( TSimul, exo_nbr );
 	exo_simul = tmp_exo_simul;
        exo_simul(:,1:3)=eps';

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

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

	Var_cons1_index = strmatch( Var_cons1, endo_names, 'exact' );
        Var_cons1_steady = muIss;
        Var_cons2_index = strmatch( Var_cons2, endo_names, 'exact' );
        Var_cons2_steady = muEss;
	SIMU_CONS = zeros( size( dr.ys, 1 ), iter+maximum_lag );
	SIMU = SIMU_CONS;
	
	oghu = zeros( size( dr.ghu ) );
	oghu( order_var, : ) = dr.ghu;
	ZIMat = zeros( 2*TshadowI, 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( 'Model7: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 i = 2:iter+maximum_lag
		yhat_cons = SIMU_CONS( order_var( k2 ), i-1 );
		yhat = SIMU( order_var( k2 ), i-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_( i-1, : )';
		epstmp = oghu * epsilon;
		
		SIMU( :, i ) = 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( M, V, AlphaStart, sDM, TshadowI, Accuracy, Options, ZeroVec, OneVec );

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

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

		
		if alpha( 1 ) > Accuracy
			CurrentBindPeriod = CurrentBindPeriod + 1;
		else
			if CurrentBindPeriod > LongestBindPeriod
				LongestBindPeriod = CurrentBindPeriod;
			end
			CurrentBindPeriod = 0;
		end
	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;
    % muihat and muehat are in row 8 and 9. Output is in row 24.
    SIMU=SIMU(:,2:end);
    SIMU_CONS=SIMU_CONS(:,2:end);
    
    
%Discard the first periods as these are affected by initial values and
%therefore may bias our results:

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

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

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

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

% Saving data in levels
SimOutputLevels(:,jj,:) = SIMU_CONS_check;

disp(jj);

end


% Saving the matrices and data
sf=sbar*100;
filename=['SimOutputSkew_sbar=' num2str(sf,2) '_pct'];
save(filename,'SimOutputLevels','sbar');

toc;
