clear;
format('compact')
global solF etah kappa A delta rho targetR currente currentK sigma m phi KKK beta lambda criticalc gamma maversion dt newequity n newe

load simulation_Baseline.mat;
%load simulationbaseline_old.mat
rng(1) ;  % It is important to have a fixed initial seed, because simulation resutls tend to vary.

e0=ec;  % Start from the cutoff of distress region.  

Requity=[-0.02
    -0.05
    -0.10
    -0.15
    -0.3];

R=Requity;
T=numel(R);e=zeros(T,1);K=zeros(T,1);investment=zeros(T,1);landprice=zeros(T,1);Sharperatio=zeros(T,1);equity=zeros(T,1);Sharpe2=zeros(T,1);
e(1)=e0;K(1)=1;
outcome=ils(e(1));
investment(1)=outcome(1)*K(1);
    landprice(1)=outcome(2)*K(1);
    Sharperatio(1)=outcome(3);
    equity(1)=min((1-lambda)*(outcome(2)+outcome(4)),e(1))*K(1);
dZ=zeros(T,1);totalshock=zeros(T,1);terme=zeros(T,1); 
n=6;
dt=0.25;t=1;length=2;xx=numel(qeprime2);num=4000;
for t=1:T
    currente=e0; currentK=1; 
    targetR=0;     initaldrift = myfunequity_multiquarters( 0 ); 
    targetR = R(t) ;
    if(t<=2)
    targetR=R(t) + initaldrift;
    end
   % targetR=R(t);
    initial_guess=0;
    dZ(t)=fsolve(@myfunequity_multiquarters,initial_guess);
    totalshock(t)=(1+sigma*dZ(t))^n-1;
    terme(t)=newe;
    simulatedLength = TT*length  ;  % Start from the end of period equity capital, and see how many periods we have crisis.
    n4sigma=0;n5sigma=0;
    
    for x=1:num
        hitting=0;awienerp=0;leng=0;k=0;
        ee=zeros(simulatedLength,1);%wi=zeros(TT*length,1);
        estartt=newe;  % Start from the end of period equity capital by the shock
        tt=1;KKstart=1;numnocrisis=0;
        % starts simulating.    
        ee(1)=estartt;indicator=ones(simulatedLength,1);KK(1)=1;equity(1)=ee(1);
        tt=2;
        while tt <= simulatedLength
            %computing mu_y and sigma_y
            yy=min(floor(1+(ee(tt-1)-elbarnew)/de),xx-1);yyy=yy+1;weight=(ee(tt-1)-(yy-1)*de-elbarnew)/de;
            mu_e=weight*mue(yyy)+(1-weight)*mue(yy);
            sigma_e=weight*sigmae(yyy)+(1-weight)*sigmae(yy);
            %Foneprime=w*Fprime(ttt)+(1-w)*Fprime(tt);
            %Fdoubleprime=w*F2prime(ttt)+(1-w)*F2prime(tt);
            dwienerp=sqrt(dt)*randn;%    awienerp=awienerp+dwienerp;
            ee(tt)=ee(tt-1)+mu_e*dt+sigma_e*dwienerp;
            ee(tt)=max(min(ee(tt),eubar),elbarnew);
            tt=tt+1;
        end
        if min(ee)<ec;%rpave*4;%0.2;
             n4sigma=n4sigma+1;
        end
        if min(ee)<criticalc;%sharpe(t-1)>=baselineSR*16;%rpave*4;%0.2;
            n5sigma=n5sigma+1;
        end
    end
    n4sigma=n4sigma/num;n5sigma=n5sigma/num;
    crisis(t)=n5sigma;Sharpe2(t)=n5sigma;%Sharpe3(j)=n3sigma;
end

totalshock
Requity./totalshock
crisis
output=[totalshock,terme,Sharpe2];
xlswrite('StressTest.xlsx', output);


terme

equities = min((1-lambda)*(pe+qe), espan)
figure; plot( espan(1:3000) ,  equities(1:3000) )

% for t=1:T
%     currente=e0;targetR=R(t); currentK=1; 
%     dZ(t)=fsolve(@myfunequity_multiquarters,0);
%     totalshock(t)=(1+sigma*dZ(t))^n-1;
% end
% totalshock = (1+sigma.*dZ).^n - 1
% Requity./totalshock
% 
% e0 = ec ;  
% for t = 1:numel(Requity)
%     currente=e0;targetR=R(t); currentK=1; 
%     %dZ(t)=fsolve(@myfunequity,0);
%     dZ(t)=fsolve(@myfunequity_multiquarters_New_2017_04,0);
% end
% totalshock = (1+sigma.*dZ).^n - 1
% Requity./totalshock



% Study the properties of myfunequity_multiquarters
e0 = ec ;    n=3; dt=0.25; currente=e0;  currentK=1;  targetRs =linspace( 0, -0.35, 30 )  ; dZs = zeros( numel(targetRs), 1 );
for i = 1:numel(targetRs)
targetR=0;     initaldrift = myfunequity_multiquarters( 0 ); 
targetR= targetRs(i) + initaldrift;
initial_guess =0;
dZs(i) = fsolve(@myfunequity_multiquarters,initial_guess);
%dZs(i) = fsolve(@myfunequity,initial_guess); 
end
totalshocks = (1+sigma.*dZs).^n - 1 ;

fig = figure;
plot( totalshocks, targetRs )
xlabel( 'Total Shock to K_t' )
ylabel(  'Change in Equity' )
%OutputPDF( fig, 'ShocktoEquity_Old.pdf' )
OutputPDF( fig, [ 'ShocktoEquity.pdf' ] )

fig = figure;
plot( abs(totalshocks) , targetRs.' ./ totalshocks )
%OutputPDF( fig, 'Amplification_Old.pdf' )
OutputPDF( fig, ['Amplification.pdf' ] )

close all;



