% simulation weekly data; annual average are calculated for statistics
format('compact')
clear;
load simulation_Baseline.mat;
rng(1) ;  % It is important to have a fixed initial seed, because simulation resutls tend to vary.

riskpremium=qe./we.*rp_k+pe./we.*rp_h;
lambdaprime=lambda + 0.022;
tthetaprime=we./espan.*(espan<criticalc)+1./(1-lambdaprime).*(espan>=criticalc);
sigmaeprime=(espan.*we.*sigma.*(m.*tthetaprime-1))./(we-espan.*m.*tthetaprime.*weprime); %checked 

musigmaeprime=espan.*(m.*r+m.*tthetaprime.*(riskpremium-r)-etah-i);
% Fe_prime=we./espan-m.*tthetaprime.*weprime; %checked
% Ge_prime=c.*kappa.*Fe_prime+gamma.*qe.*qeprime.*m.*(1-tthetaprime).*we; % checked
% musigmaeprime=c.*kappa./Ge_prime.*((m.*(1-ttheta).*(rho+gamma.*i-gamma.*(0.5.*(qe.*qeprime2+qeprime.^2).*sigmae.^2)./(c.*kappa)...
%     -0.5*gamma.*(1+gamma).*(c.*sigma-qe.*qeprime.*sigmaeprime./kappa).^2./c.^2)-i-eta).*we...
%     +m.*tthetaprime.*(0.5.*(weprime2).*sigmae.^2+(A-delta)/phi+(pe-(1-phi)/phi).*i-(1-phi)/phi.*kappa.*0.5.*(i.^2)+delta.*(1-qe)));  % rewritten, checked
mueprime=musigmaeprime-sigma.*sigmaeprime;

TT=4;dt=1/TT;
num=4000;
length1=500;llength=[1;2;5];KKlbar=0.01; length=llength(3);
K=500;counter=0;position=zeros(K,1);positionnew=zeros(K,1);deltae=(eubar-elbarnew)/K;
checkrp=0;
num_con=0;n1sigma=0;n2sigma=0;n3sigma=0;n4sigma_1=0;n5sigma_1=0;n4sigma_2=0;n5sigma_2=0;n4sigma_3=0;n5sigma_3=0;
p_time=0;p=0;d=0;averleng=0;

xx=numel(qeprime2);num_con=0;num_un=0;xxx=0;
estartt=ec;

tic;
for x=1:num
t=1;KKstart=1;numnocrisis=0;
%display(  [  '. Progress: ' , num2str(100*x/num), '%. Time used: ',  num2str(toc/60,2), ' minutes'   ]   )
% starts simulating.    
ee(1)=estartt;indicator=ones(TT*length,1);KK(1)=1;equity(1)=ee(1);
t=2;
    while t<TT*length+1
    tt=min(floor(1+(ee(t-1)-elbarnew)/de),xx-1);ttt=tt+1;weight=(ee(t-1)-(tt-1)*de-elbarnew)/de;
    mu_e=weight*mueprime(ttt)+(1-weight)*mueprime(tt);
    sigma_e=weight*sigmaeprime(ttt)+(1-weight)*sigmaeprime(tt);
    dwienerp=sqrt(dt)*randn;%    
    ee(t)=ee(t-1)+mu_e*dt+sigma_e*dwienerp;
    ee(t)=max(min(ee(t),eubar),elbarnew);

    t=t+1;
    end
if min(ee(1:TT*llength(1)))<ec;%rpave*4;%0.2;
         n4sigma_1=n4sigma_1+1;
     end

if min(ee(1:TT*llength(1)))<criticalc;%sharpe(t-1)>=baselineSR*16;%rpave*4;%0.2;
        
        n5sigma_1=n5sigma_1+1;
     end
    
if min(ee(1:TT*llength(2)))<ec;%rpave*4;%0.2;
         n4sigma_2=n4sigma_2+1;
     end

if min(ee(1:TT*llength(2)))<criticalc;%sharpe(t-1)>=baselineSR*16;%rpave*4;%0.2;
        
        n5sigma_2=n5sigma_2+1;
end

if min(ee(1:TT*llength(3)))<ec;%rpave*4;%0.2;
         n4sigma_3=n4sigma_3+1;
     end

if min(ee(1:TT*llength(3)))<criticalc;%sharpe(t-1)>=baselineSR*16;%rpave*4;%0.2;
        
        n5sigma_3=n5sigma_3+1;
end
end

n4sigma_1=(n4sigma_1-1)/num;n5sigma_1=(n5sigma_1-1)/num;
n4sigma_2=(n4sigma_2-1)/num;n5sigma_2=(n5sigma_2-1)/num;
n4sigma_3=(n4sigma_3-1)/num;n5sigma_3=(n5sigma_3-1)/num;

crisis(1)=n5sigma_1;crisis(2)=n5sigma_2;crisis(3)=n5sigma_3;
Sharpe2(1)=n4sigma_1;Sharpe2(2)=n4sigma_2;Sharpe2(3)=n4sigma_3;
crisis  % Probability of crisis over 1, 2, and 5 years


