% simulation weekly data; annual average are calculated for statistics
clear;
load simulation_Baseline.mat;
TT=4;dt=1/TT;
num=4000;   % Number of simulations in each year, and each starting e.
length1=500;
llength=[2;5;10];   % Length of 2, 5, and 10 years
KKlbar=0.01;
etest=[criticalc-0.05:0.1:4]';   % Starting values of e

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

% simulate 500 times for each starting value of e
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=0;n5sigma=0;p_time=0;p=0;d=0;averleng=0;

% real activities
e_cgrowth=0;e_cvol=0; % consumption
e_investment=0;e_investmentvol=0; %investment
e_cgrowth_con=0;e_cgrowth_un=0;e_cvol_con=0;e_cvol_un=0;
e_investment_un=0;e_investment_con=0;e_investmentvol_un=0;e_investmentvol_con=0;
e_invvol=0;e_invvol_un=0;e_invvol_con=0;
e_pvol=0;e_pvol_un=0;e_pvol_con=0;
e_rentalpriceratio=0;e_rentalpriceratio_un=0;e_rentalpriceratio_con=0;

%asset prices
e_alphak=0;e_alphak_con=0;e_alphak_un=0;e_alphah=0;e_alphah_con=0;e_alphah_un=0;

e_rpk=0;e_rpk_con=0;e_rpk_un=0;e_rph=0;e_rph_con=0;e_rph_un=0;
std_serk=0;std_serk_un=0;std_serk_con=0;
std_serh=0;std_serh_un=0;std_serh_con=0;
std_qq=0;std_qq_con=0;std_qq_un=0;
e_rf=0;e_rf_con=0;e_rf_un=0;

e_sigmak=0;e_sigmak_un=0;e_sigmak_con=0;e_sigmah=0;e_sigmah_un=0;e_sigmah_con=0;
e_sharpe=0;e_sharpe_un=0;e_sharpe_con=0;
e_ppoverqq=0;e_qq=0;e_qq_con=0;e_qq_un=0;
%constraint
e_theta=0;std_theta=0;e_theta_un=0;e_theta_con=0;std_theta_un=0;std_theta_con=0;
e_ee=0;e_median=0;std_ee=0;e_ee_con=0;e_ee_un=0;std_ee_con=0;std_ee_un=0;
cov_equity_equity_con=0;cov_equity_equity_un=0;
cov_equity_inv_con=0;cov_equity_inv_un=0;
cov_equity_sharpe_con=0;cov_equity_sharpe_un=0;
cov_equity_consumption_con=0;cov_equity_consumption_un=0;
cov_equity_landprice_con=0;cov_equity_landprice_un=0;
cov_sharpe_sharpe_con=0;cov_sharpe_sharpe_un=0;
cov_sharpe_inv_con=0;cov_sharpe_inv_un=0;
cov_sharpe_consumption_con=0;cov_sharpe_consumption_un=0;
cov_sharpe_landprice_con=0;cov_sharpe_landprice_un=0;
cov_inv_inv_con=0;cov_inv_inv_un=0;
cov_inv_landprice_con=0;cov_inv_landprice_un=0;
cov_inv_consumption_con=0;cov_inv_consumption_un=0;
cov_consumption_consumption_con=0;cov_consumption_consumption_un=0;
cov_consumption_landprice_con=0;cov_consumption_landprice_un=0;
cov_landprice_landprice_con=0;cov_landprice_landprice_un=0;

mean_equity_growth_con=0;mean_inv_growth_con=0;mean_sharpe_con=0;mean_consumption_growth_con=0;mean_landprice_growth_con=0;
mean_equity_growth_un=0;mean_inv_growth_un=0;mean_sharpe_un=0;mean_consumption_growth_un=0;mean_landprice_growth_un=0;

xx=numel(qeprime2);num_con=0;num_un=0;xxx=0;
J=numel(etest);   
estartt=zeros(J,1);  

tic;
for m=1:3      % For cases of 2,5, and 10 years
    length=llength(m);
    for j=1:J        % For different starting values of e
        display(  [  '. Progress: ' , num2str(100*(j+(m-1)*J)/(3*J)), '%. Time used: ',  num2str(toc/60,2), ' minutes'   ]   )
        for x=1:num    % For different numbers of simulations
            hitting=0;awienerp=0;leng=0;k=0;
            ee=zeros(TT*length,1);%wi=zeros(TT*length,1);

            cgrowth=zeros(TT*length,1);cvol=zeros(TT*length,1);KK=zeros(TT*length,1);
            investment=zeros(TT*length,1);invest=zeros(TT*length,1);investmentvol=zeros(TT*length,1);inv_growth=zeros(TT*length,1);
            netinvestment=zeros(TT*length,1);
            pricevol=zeros(TT*length,1);
            rentalpriceratio=zeros(TT*length,1);indi=zeros(TT*length,1);equity=zeros(TT*length,1);equity_growth=zeros(TT*length,1);
            landprice_growth=zeros(TT*length,1);
            consumption=zeros(TT*length,1);consumption_growth=zeros(TT*length,1);

            alphak=zeros(TT*length,1);alphah=zeros(TT*length,1);
            rpk=zeros(TT*length,1);rph=zeros(TT*length,1);
            rf=zeros(TT*length,1);
            sigmak=zeros(TT*length,1);sigmah=zeros(TT*length,1);
            sharpe=zeros(TT*length,1);

            theta=zeros(TT*length,1);

            serk=zeros(TT*length,1);serh=zeros(TT*length,1);
            ppoverqq=zeros(TT*length,1);
            qq=zeros(TT*length,1);pp=zeros(TT*length,1);

            estartt(j)=etest(j);

            t=1;KKstart=1;numnocrisis=0;

            % starts simulating.    
            ee(1)=estartt(j);indicator=ones(TT*length,1);KK(1)=1;equity(1)=ee(1);
            t=2;
            % In each simulation, we have to calculate ee for TT*length+1 quarters, where length is number of years, TT is number of time points in each year
            while t<TT*length+1   
            %computing mu_y and sigma_y
            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*mue(ttt)+(1-weight)*mue(tt);
            sigma_e=weight*sigmae(ttt)+(1-weight)*sigmae(tt);
            %Foneprime=w*Fprime(ttt)+(1-w)*Fprime(tt);
            %Fdoubleprime=w*F2prime(ttt)+(1-w)*F2prime(tt);
            cgrowth(t-1)=weight*muc(ttt)+(1-weight)*muc(tt);cvol(t-1)=weight*sigmac(ttt)+(1-weight)*sigmac(tt);
            invest(t-1)=(weight*inv(ttt)+(1-weight)*inv(tt))*KK(t-1);
            investment(t-1)=(weight*inv(ttt)+(1-weight)*inv(tt));
            netinvestment(t-1)=(weight*(i(ttt)+delta)+(1-weight)*(i(tt)+delta));
            investmentvol(t-1)=(weight*Ivol(ttt)+(1-weight)*Ivol(tt))+sigma;
            pricevol(t-1)=weight*pvol(ttt)+(1-weight)*pvol(tt);
            rentalpriceratio(t-1)=weight*Doverp(ttt)+(1-weight)*Doverp(tt);

            alphak(t-1)=weight*alpha_k(ttt)+(1-weight)*alpha_k(tt);alphah(t-1)=weight*alpha_h(ttt)+(1-weight)*alpha_h(tt);
            qq(t-1)=weight*qe(ttt)+(1-weight)*qe(tt);pp(t-1)=weight*pe(ttt)+(1-weight)*pe(tt);
            sigmak(t-1)=weight*sigmaRk(ttt)+(1-weight)*sigmaRk(tt);sigmah(t-1)=weight*sigmaRh(ttt)+(1-weight)*sigmaRh(tt);
            rpk(t-1)=weight*rp_k(ttt)+(1-weight)*rp_k(tt);rph(t-1)=weight*rp_h(ttt)+(1-weight)*rp_h(tt);
            sharpe(t-1)=weight*SR(ttt)+(1-weight)*SR(tt); %sharpe ratio
            theta(t-1)=weight*ttheta(ttt)+(1-weight)*ttheta(tt);
            rf(t-1)=weight*r(ttt)+(1-weight)*r(tt);
            ppoverqq(t-1)=pp(t-1)/(pp(t-1)+qq(t-1));
            consum_scaled(t-1)=weight*c(ttt)+(1-weight)*c(tt);
            consumption(t-1)=consum_scaled(t-1)*KK(t-1); %consumption levels
            landprice(t-1)=pp(t-1)*KK(t-1); %land prices

            equity(t-1)=min(ee(t-1),(1-lambda)*(pp(t-1)+qq(t-1)))*KK(t-1); %bank equity
                srave=0.03;    
            baselineSR=0.3832226936;
            if sharpe(t-1)>=baselineSR*1;
                n1sigma=n1sigma+1;
            end
            if sharpe(t-1)>=baselineSR*2; 
            n2sigma=n2sigma+1;
            end
            if sharpe(t-1)>=baselineSR*4;%0.1;
                n3sigma=n3sigma+1;
            end
            dwienerp=sqrt(dt)*randn;%    awienerp=awienerp+dwienerp;
            ee(t)=ee(t-1)+mu_e*dt+sigma_e*dwienerp;
            ee(t)=max(min(ee(t),eubar),elbarnew);

            KK(t)=KK(t-1)*(1+netinvestment(t-1)*dt-delta*dt+sigma*dwienerp);
            KK(t)=max(KK(t),KKlbar);

            serk(t-1)=rpk(t-1)*dt+sigmak(t-1)*dwienerp;
            serh(t-1)=rph(t-1)*dt+sigmah(t-1)*dwienerp;

            t=t+1;
            end
            if min(ee)<ec;     % To calculate whether ee is below ec, i.e. whether there is a distress
                    n4sigma=n4sigma+1;
            end
            if min(ee)<criticalc;   % To calculate whether ee is below criticalc, i.e. whether there is a crisis.
                    n5sigma=n5sigma+1;
            end
        end
        n1sigma=n1sigma/num/TT/length;n2sigma=n2sigma/num/TT/length;
        n3sigma=n3sigma/num/TT/length;n4sigma=(n4sigma-1)/num;n5sigma=(n5sigma-1)/num;
        distress(j,m)=n5sigma;Sharpe2(j,m)=n4sigma;
    end
end

save simulation_Baseline_ProbCrisis_full.mat
save('simulation_Baseline_ProbCrisis.mat',  'distress', 'estartt');

Figure9 = figure;
subplot(121)
plot(estartt,Sharpe2(:,1), '-b'   , estartt,Sharpe2(:,2), '-.b', estartt,Sharpe2(:,3), '--b');
title('Probability of being distressed')
legend('in next 2 years',1,'in next 5 years',2,'in next 10 years',3)
xlabel('starting value e_{init}')
subplot(122)
plot(estartt,distress(:,1), '-b'   , estartt,distress(:,2), '-.b', estartt,distress(:,3), '--b');
title('Probability of capital constraint being binding')
legend('in next 2 years',1,'in next 5 years',2,'in next 10 years',3)
xlabel('starting value e_{init}')
cd('/Users/wenhao/Dropbox/housing/MatlabCodes_and_Figures_ModifiedbyWenhao_2017_02/Figures and Tables for Paper')
OutputPDF( Figure9, 'Figure9', [8,4] );



