% Henrik Jensen, Søren Hove Ravn, Emiliano Santoro
% Changing Credit Limits, Changing Business Cycles

% Function, used to solve for "shadow prices"
% Based on the work of Tom Holden's paper: http://www.dynare.org/wp-repo/dynarewp004.pdf


TShadow=20;
TIRF=20;

oo_.steady_state(1:33,1)=steadys;

Options = optimset( 'MaxIter', 4096, 'LargeScale', 'on' );
ZeroVec = zeros( 2*TShadow, 1 );

shock1=exist('ahat_epsilon_eahat');
shock2=exist('ahat_epsilon_eshat');
shock3=exist('ahat_epsilon_epshat');

if shock1==1;
    shock = 'eahat';
end
if shock2==1;
    shock = 'eshat';
end
if shock3==1;
    shock = 'epshat';
end

		alpha = ZeroVec;


			gIe = GetIRFnew( 'muehat', shock, TIRF );
            gIi = GetIRFnew( 'muihat', shock, TIRF );
            gI = [gIi(1:TShadow) ; gIe(1:TShadow)];

			for indx = 0:( TShadow - 1 )
				gIShadowee( :, indx + 1 ) = muEss.*GetRelIRF( 'muehat', [ 'shadowEnt_' sprintf( '%d', indx ) ]);
                gIShadowei( :, indx + 1 ) = muEss.*GetRelIRF( 'muehat', [ 'shadowImp_' sprintf( '%d', indx ) ]);
                gIShadowii( :, indx + 1 ) = muIss.*GetRelIRF( 'muihat', [ 'shadowImp_' sprintf( '%d', indx ) ]);
                gIShadowie( :, indx + 1 ) = muIss.*GetRelIRF( 'muihat', [ 'shadowEnt_' sprintf( '%d', indx ) ]);
            end
            
            gIShadow = [ gIShadowii(1:TShadow,:) gIShadowie(1:TShadow,:) ; gIShadowei(1:TShadow,:) gIShadowee(1:TShadow,:) ];

            % Now solve the optimization problem as described in the paper
            % by Holden & Paetz. The problem is:
			% alpha' * ( gI( 1:TShadow ) + gIShadow( 1:TShadow, : ) * alpha ) = 0
			% => 1/2 * alpha' * ( gIShadow( 1:TShadow, : ) + gIShadow( 1:TShadow, : )' ) * alpha + gI' * alpha = 0
			% gI + gIShadow * alpha >= 0
			% => -gIShadow * alpha <= gI
            
            % This problem can be solved using the quadprog function. See
            % Matlab Help for details of the inputs etc.

			[ alpha, FoundValue ] =  quadprog( gIShadow + gIShadow', gI, -gIShadow, gI, [], [], ZeroVec, [], alpha, Options );

			% gI = gI + gIShadow * alpha;
            
            % Finally, check that the value of the function is in fact
            % zero, then use the argument to compute and plot bounded irfs

			if FoundValue < 1e-5
                fprintf( 1, [ 'Solved problem with relevant shock; Value=' ], sum( alpha > 1e-15 ) );
                fprintf('%d\n',FoundValue);

            y_bound = GetBoundedIRF_bothnew('yhat',shock,alpha,TIRF);
            i_bound = GetBoundedIRF_bothnew('ihat',shock,alpha,TIRF);
            mui_bound = GetBoundedIRF_bothnew('muihat',shock,alpha,TIRF);
            mue_bound = GetBoundedIRF_bothnew('muehat',shock,alpha,TIRF);
            qh_bound = GetBoundedIRF_bothnew('qhat',shock,alpha,TIRF);
            qk_bound = GetBoundedIRF_bothnew('qkhat',shock,alpha,TIRF);
            r_bound = GetBoundedIRF_bothnew('rhat',shock,alpha,TIRF);
            cp_bound = GetBoundedIRF_bothnew('cphat',shock,alpha,TIRF);
            ci_bound = GetBoundedIRF_bothnew('cihat',shock,alpha,TIRF);
            ce_bound = GetBoundedIRF_bothnew('cehat',shock,alpha,TIRF);
            hp_bound = GetBoundedIRF_bothnew('hphat',shock,alpha,TIRF);
            hi_bound = GetBoundedIRF_bothnew('hihat',shock,alpha,TIRF);
            he_bound = GetBoundedIRF_bothnew('hehat',shock,alpha,TIRF);
            bp_bound = GetBoundedIRF_bothnew('bphat',shock,alpha,TIRF);
            bi_bound = GetBoundedIRF_bothnew('bihat',shock,alpha,TIRF);
            be_bound = GetBoundedIRF_bothnew('behat',shock,alpha,TIRF);
            np_bound = GetBoundedIRF_bothnew('nphat',shock,alpha,TIRF);
            ni_bound = GetBoundedIRF_bothnew('nihat',shock,alpha,TIRF);
            wp_bound = GetBoundedIRF_bothnew('wphat',shock,alpha,TIRF);
            wi_bound = GetBoundedIRF_bothnew('wihat',shock,alpha,TIRF);
            ctot_bound = GetBoundedIRF_bothnew('ctot',shock,alpha,TIRF);
            ntot_bound = GetBoundedIRF_bothnew('ntot',shock,alpha,TIRF);


            y_unbound = GetIRFnew('yhat',shock,TIRF);
            i_unbound = GetIRFnew('ihat',shock,TIRF);
            mui_unbound = GetIRFnew('muihat',shock,TIRF);
            mue_unbound = GetIRFnew('muehat',shock,TIRF);
            qh_unbound = GetIRFnew('qhat',shock,TIRF);
            qk_unbound = GetIRFnew('qkhat',shock,TIRF);
            r_unbound = GetIRFnew('rhat',shock,TIRF);
            cp_unbound = GetIRFnew('cphat',shock,TIRF);
            ci_unbound = GetIRFnew('cihat',shock,TIRF);
            ce_unbound = GetIRFnew('cehat',shock,TIRF);
            hp_unbound = GetIRFnew('hphat',shock,TIRF);
            hi_unbound = GetIRFnew('hihat',shock,TIRF);
            he_unbound = GetIRFnew('hehat',shock,TIRF);
            bp_unbound = GetIRFnew('bphat',shock,TIRF);
            bi_unbound = GetIRFnew('bihat',shock,TIRF);
            be_unbound = GetIRFnew('behat',shock,TIRF);
            np_unbound = GetIRFnew('nphat',shock,TIRF);
            ni_unbound = GetIRFnew('nihat',shock,TIRF);
            wp_unbound = GetIRFnew('wphat',shock,TIRF);
            ctot_unbound = GetIRFnew('ctot',shock,TIRF);
            ntot_unbound = GetIRFnew('ntot',shock,TIRF);

			else
				fprintf( 1, [ 'Found minimum that is not exactly zero; Value=' ], FoundValue );
                fprintf('%d\n',FoundValue);
            end
		    
            xx=(0:TIRF-1)';
%             
%             figure;
%             subplot(2,2,1)
%             plot(xx,y_unbound(1:TIRF),'-','MarkerSize',7,'LineWidth',1.3);
%             hold on;
%             plot(xx,GetMean('yhat')*ones(TIRF),'-k','MarkerSize',3,'LineWidth',1.3);
%             plot(xx,y_bound(1:TIRF),'r','MarkerSize',7,'LineWidth',1.3);
%             title('Output Level','fontsize',11);
%             hold off;
%             subplot(2,2,2)
%             plot(xx,i_unbound(1:TIRF),'-','MarkerSize',7,'LineWidth',1.3);
%             hold on;
%             plot(xx,GetMean('ihat')*ones(TIRF),'-k','MarkerSize',3,'LineWidth',1.3);
%             plot(xx,i_bound(1:TIRF),'r','MarkerSize',7,'LineWidth',1.3);
%             title('Investment Level','fontsize',11);  
%             hold off;
%             subplot(2,2,3)
%             plot(xx,mui_unbound(1:TIRF),'-','MarkerSize',7,'LineWidth',1.3);
%             hold on;
%             plot(xx,GetMean('muihat')*ones(TIRF),'-k','MarkerSize',3,'LineWidth',1.3);
%             plot(xx,mui_bound(1:TIRF),'r','MarkerSize',7,'LineWidth',1.3);
%             title('Shadow Price I Level','fontsize',11);  
%             hold off;
%             subplot(2,2,4)
%             plot(xx,mue_unbound(1:TIRF),'-','MarkerSize',7,'LineWidth',1.3);
%             hold on;
%             plot(xx,GetMean('muehat')*ones(TIRF),'-k','MarkerSize',3,'LineWidth',1.3);
%             plot(xx,mue_bound(1:TIRF),'r','MarkerSize',7,'LineWidth',1.3);
%             title('Shadow Price E Level','fontsize',11);  

            AA=[ (y_bound-Yss)./Yss (i_bound-Iss)./Iss (mui_bound-muIss)./muIss (mue_bound-muEss)./muEss (qh_bound-Qss)./Qss (r_bound-Rss)./Rss ...
                (cp_bound-CPss)./CPss (ci_bound-CIss)./CIss (ce_bound-CEss)./CEss (hp_bound-HPss)./HPss (hi_bound-HIss)./HIss (he_bound-HEss)./HEss ...
                (bp_bound-BPss)./BPss (bi_bound-BIss)./BIss (be_bound-BEss)./BEss (np_bound-NPss)./NPss (ni_bound-NIss)./NIss (wp_bound-WPss)./WPss (wi_bound-WIss)./WIss];
            
            AA_report= [ (y_bound-Yss)./Yss (ctot_bound-Ctotss)./Ctotss (i_bound-Iss)./Iss -((bp_bound-BPss)./BPss) (qh_bound-Qss)./Qss (r_bound-Rss)./Rss  (ntot_bound-Ntotss)./Ntotss];
            AA_report_large = [(y_bound-Yss)./Yss (y_unbound-Yss)./Yss];
            
            BA=(y_unbound-Yss);
            BB=BA./Yss;
            nulline=zeros(1,TIRF);
%             figure
%             plot(xx,AA(1:TIRF,1),'-','MarkerSize',7,'LineWidth',1.3);
%             hold on;
%             plot(xx,nulline,'-k','MarkerSize',3,'LineWidth',1.3);
%             plot(xx,BB(1:TIRF),'r','MarkerSize',7,'LineWidth',1.3);
%             title('Output IRF','fontsize',11);
%             hold off;
          

