%% ********** Section: Initialization ********** 
format('compact')
clear
%% ********** Section: Initialization ********** 
global rho sigma sigmasqr lambda gamma A phi r elbarnew wlbarnew qinf pinf plbarnew...
    qlbarnew eubar kappa delta etah etal eentry m iota tau uppersharpe maversion solF ...
    zeta chi  elbarnewH
% Key Paramters for Matching Moments
etah=0.12;  etal=etah;  % We have eta = etah = etal
lambda=0.75; % ===== This is different from the baseline. =====
A=0.133;  % productivity.  Target the Investment-to-capital ratio to be 11%
phi=0.4;    
gamma=0.13;   %gamma is the effective risk aversion in the Euler equation for the interest rate. Equation (1) in Derivation072016. Target interest rate volatility of 1%
tau=  0.88;  % This tau is just tau=p(elbar)beta/(1+elbar*beta). Old value: 0.419. 
uppersharpe=6.5;  % Entry barrier parameter B
maversion=1.8;

% Other parameters
kesi=tau;y=0;
 rho=0.02;
 sigma=0.03;sigmasqr=sigma^2;
kappa=3;delta=0.1;
 m=1;  % m is abandoned in the new version of the notes. Instead, maversion is used. It should be set to 1.
iota=0;
 eentry=4;
 zeta = 1;   % !!! Newly added by Wenhao
 chi = 0.85;  % !!! Newly added by Wenhao. The value of chi actually doesn't matter at all, as a property of the model.
qinf=(-(rho+delta+(((2*maversion-gamma*(1+gamma)*(1-lambda))*sigmasqr/(2*(1-lambda))))-(gamma/kappa))...
    +sqrt((rho+delta+(((2*maversion-gamma*(1+gamma)*(1-lambda))*sigmasqr/(2*(1-lambda))))-(gamma/kappa))^2+4*gamma*A/kappa))/(2*gamma/kappa);
iinf=((qinf-1)/kappa);   
iinfgross=delta + iinf+0.5*kappa*iinf^2;
pinf=(((1-phi)/phi)*(A-iinfgross))/(rho+(gamma-1)*iinf+((2*maversion-gamma*(1+gamma)*(1-lambda))*sigmasqr/(2*(1-lambda))));
rinf=rho+gamma*iinf-0.5*gamma*(1+gamma)*sigmasqr; %CF=A-kappa/2*iinf^2;  See the first formula in A.5.1 in the old version of the paper
winf=qinf+pinf;eubar=20;
elbarnewH = winf/( uppersharpe*(1-m*tau)/sigma/maversion + tau  );  % This is the upper bound of elbarnew that makes winf >= wlbarnew
dif=1;%dif=winf-wlbarnew;
KK=500;


%% **********  Section: Initial Solution ********** 
% General idea: Evaluate the ODE program over a grid of elbarnew and plbarnew. Then narrow the range of the grid. Finally solve the proble in a very niche grid.
% Note the vector is y(x) = [ p(x), q(x), p'(x), q'(x)]. The technique of continuation is used.
options=bvpset('relTol',1e-3,'NMax',300000,'Vectorized','on');%,'BCJacobian',@bcjacnew);
options2=bvpset('relTol',1e-12,'NMax',300000,'Vectorized','on');%,'BCJacobian',@bcjacnew);  % Higher accuracy. Default is 1e-3

% ***** Evaluate certain values and initialize the whole problem.  ***** 
% Replicate the old results
elbarnew = 0.2; qlbarnew = 0.9 ; wlbarnew=uppersharpe*elbarnew*(1-m*tau)/sigma/maversion+elbarnew*tau;
plbarnew = wlbarnew - qlbarnew; beta=tau/(plbarnew-tau*elbarnew) ;  
solinitF= bvpinit(linspace(elbarnew,eubar,KK),@Finitwithproduction);  
solF = bvp4c(@Fvecwithproduction,@Fbcwithproduction,solinitF,options);
solinitF= bvpinit(linspace(elbarnew,eubar,KK),@Finitwithproduction2);  
solF=bvp4c(@Fvecwithproduction,@Fbcwithproduction,solinitF,options2);
mis=deval(solF,elbarnew,3)-tau-deval(solF,elbarnew,4)            
subplot(2,1,1);   plot( solF.x, solF.y(1,:) )
subplot(2,1,2);   plot( solF.x, solF.y(2,:) )

%% Search qlbarnew as a region
% %*********************** First round: Find the rough region of possible solution ***********************
N = 10;  Diffs = nan*ones(N,N); 
elbarnews = linspace(0.04, elbarnewH, N);
qlbarnews = linspace(0.75, 1, N) ;

while(0)
    for i = 1:N
        elbarnew = elbarnews(i);
        display( ['******  New Round:  elbar = ', num2str(elbarnew) , '  *******'] )
        wlbarnew=uppersharpe*elbarnew*(1-m*tau)/sigma/maversion+elbarnew*tau;
        for j = 1:N
           display(j)
           qlbarnew=qlbarnews(j);
           plbarnew=wlbarnew-qlbarnew;
           if( plbarnew < pinf && qlbarnew < qinf && plbarnew>0 && qlbarnew>0  && plbarnew>tau*elbarnew )  % Only calculate cases with increasing q(e) and p(e)    
               try
                    solinitF= bvpinit(linspace(elbarnew,eubar,KK),@Finitwithproduction2);  
                    solF=bvp4c(@Fvecwithproduction,@Fbcwithproduction,solinitF,options);
                    mis=deval(solF,elbarnew,3)-tau-deval(solF,elbarnew,4)  % deval(solF,elbarnew,3) returns the 3rd component of the evaluated solution
                    Diffs(i,j) = mis ;
               catch ME
                    %display(ME)
               end 
           end
        end
    end
    OutputSimulation(elbarnews, qlbarnews, Diffs)
    load chirp; sound(y,Fs)   % Produce a sound to signal the end of the loop
% ------ End Loop  ------ 
end

%*** More rounds: Find A Finer region of possible solution. Iterate these again and again. ***
N = 10;  Diffs = nan*ones(N,N); 
elbarnews = linspace(0.087, 0.089, N);
qlbarnews = linspace(0.85, 1, N) ;  % note: adjust rnages so that err will alternative between negative and positive values.
    
% ****************************** Bisection Method to Solve the Problem ******************************
mistol=1e-3;errtol=1e-3; earlytol = 1e-5;
elbar0 = min(elbarnews);   elbar00 = max(elbarnews) ;      % !!! Added by Wenhao, to enlarge the range of search.
elbarnew= (elbar0+elbar00)/2;
errors = []; elbars = []; qlbars = [];

% Step 1: Initialize the err variable. The initialization is very very important because it determines the most important first partition. If it is wrong, then there is no possibility of finding an accurate solution.
mis=100;  mis_prev = mis+1;
qlbar0=qlbarnews(1); qlbar00=qlbarnews(N);qlbarnew=(qlbar0+qlbar00)/2;  % Initialize qlbar in each round of elbar
display(['********** New elbar = ', num2str(elbarnew), ' *********** ' ])
while abs(mis)>mistol  && abs(mis - mis_prev) > mistol
        plbarnew=wlbarnew-qlbarnew;
        solinitF= bvpinit(linspace(elbarnew,eubar,KK),@Finitwithproduction2);
        solF=bvp4c(@Fvecwithproduction,@Fbcwithproduction,solinitF,options2);    
        mis_prev = mis;
        mis=deval(solF,elbarnew,3)-tau-deval(solF,elbarnew,4);
        display([ 'qlbarnew=', num2str(qlbarnew)  ,'.  mis = ', num2str(mis)])
        if mis<0
            qlbar0=qlbarnew;qlbar00=qlbar00;qlbarnew=(qlbar0+qlbar00)/2;
        else
            qlbar0=qlbar0;qlbar00=qlbarnew;qlbarnew=(qlbar0+qlbar00)/2;
        end
end
err=deval(solF,elbarnew,4)
err_prev = err + 10*errtol;  % Initialize previous error, so that the difference is larger than errtol.
% Step 2: Main Round
while abs(err)>errtol && abs(err_prev - err) > earlytol  % Outer loop: search elbar until q'(elbar) is close to zero enough
    if err>0
        elbar0=elbarnew;elbar00=elbar00;elbarnew=(elbar0+elbar00)/2;
    else
        elbar0=elbar0;elbar00=elbarnew;elbarnew=(elbar0+elbar00)/2;
    end 
    display(['********** New elbar = ', num2str(elbarnew), ' *********** ' ])
    wlbarnew=uppersharpe*elbarnew*(1-m*tau)/sigma/maversion+elbarnew*tau;    
    dif=winf-wlbarnew;
    qlbar0=qlbarnews(1); qlbar00= min(wlbarnew, qlbarnews(N)) ;qlbarnew=(qlbar0+qlbar00)/2;  % Initialize qlbar in each round of elbar
    plbarnew=wlbarnew-qlbarnew;
    solinitF= bvpinit(linspace(elbarnew,eubar,KK),@Finitwithproduction2);
    solF=bvp4c(@Fvecwithproduction,@Fbcwithproduction,solinitF,options2);
    mis=deval(solF,elbarnew,3)-tau-deval(solF,elbarnew,4) ;
    mis_prev = mis + 10*mistol;  % Keep record of the previous mis so that I can judge whether the loop is making progress.
    while abs(mis)>mistol && abs(mis_prev - mis) > mistol
        if mis<0
            qlbar0=qlbarnew;qlbar00=qlbar00;qlbarnew=(qlbar0+qlbar00)/2;
        else
            qlbar0=qlbar0;qlbar00=qlbarnew;qlbarnew=(qlbar0+qlbar00)/2;
        end
        plbarnew=wlbarnew-qlbarnew;
        solinitF= bvpinit(linspace(elbarnew,eubar,KK),@Finitwithproduction2);
        try
            solF=bvp4c(@Fvecwithproduction,@Fbcwithproduction,solinitF,options2);    
            % If no error, continue evaluating the codes below
            mis_prev = mis;
            mis=deval(solF,elbarnew,3)-tau-deval(solF,elbarnew,4);
            display([ 'qlbarnew=', num2str(qlbarnew)  ,'.  mis = ', num2str(mis)])
        catch
        end
    end
     err_prev = err;
     err=deval(solF,elbarnew,4);
     % Question: what is the use of y ? 
     display(['Miniumum error at elbar=', num2str(elbarnew), ' is ',  num2str(err)])
     errors = [errors; err];
     elbars = [ elbars; elbarnew  ]; 
     qlbars = [ qlbars; qlbarnew ];
end
mis_to_store = mis;  err_to_store = err;
scatter(  elbars, errors )
fitvars = polyfit(elbars, errors, 1);  elbar_opt = - fitvars(2) / fitvars(1) ;    % See how errors change, and fit a line to get the intersection point with elbar axes.
display(['Next guess of elbarnew is ',  num2str(elbar_opt)])
load chirp; sound(y,Fs); 


%% **********  Section: Calculate Key Variables and Output  **********
for(i=1)  % This for is only to wrap codes
K=20000;
eupper=eubar;
de=(eupper-elbarnew)/K;
espan=[elbarnew:de:eupper];
%espan1=[elbarnew+de:de:eupper];
%espan=solF.x;
pe=deval(solF,espan,1);
beta=tau/(pe(1)-tau*elbarnew) ;   % Question: What is the correct beta?
[prime1,prime2]=deval(solF,espan,3);
peprime=prime1;
peprime2=prime2;
qe=deval(solF,espan,2);
[prime1,prime2]=deval(solF,espan,4);
qeprime=prime1;
qeprime2=prime2;
we=zeta*pe+qe;weprime=zeta*peprime+qeprime;weprime2=zeta*peprime2+qeprime2;
eta=etal.*(espan<eentry)+etah.*(espan>=eentry);
i=(qe-1)/kappa;
inv=i+delta;%+((qe-1)/kappa).^2/2*kappa;
c=A-delta-i-i.^2/2*kappa;
ttheta=max(we./espan,1/(1-lambda)); %checked
alpha_k=ttheta.*qe./we;alpha_h=ttheta.*pe./we;
sigmae=(espan.*we.*sigma.*(m.*ttheta-1))./(we-espan.*m.*ttheta.*weprime); %checked 

Ivol=qeprime.*sigmae/kappa./inv;
pvol=peprime.*sigmae./pe;
Fe=we./espan-m.*ttheta.*weprime; %checked
Ge=c.*kappa.*Fe+gamma.*qe.*qeprime.*m.*(1-ttheta).*we; % checked
musigmae=c.*kappa./Ge.*((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.*sigmae./kappa).^2./c.^2)-i-eta).*we...
    +m.*ttheta.*(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

mue=musigmae-sigma.*sigmae;
mueper=(musigmae-sigma.*sigmae)./espan; % this is correct
rr=(rho+gamma.*i).*c.*kappa.*(we./espan-m.*ttheta.*(weprime))...
    -gamma.*qe.*qeprime.*((-i-eta).*we+m.*ttheta.*(0.5.*weprime2.*sigmae.^2+(A-delta)./phi+(pe-((1-phi)./phi)).*i-((1-phi)./phi).*0.5.*kappa.*i.^2+delta.*(1-qe)))...
    -gamma.*(qe.*qeprime2+qeprime.^2).*sigmae.^2.*0.5.*(we./espan-m.*ttheta.*(weprime))...
    -0.5*gamma.*(1+gamma).*(we./espan-m.*ttheta.*(weprime)).*kappa.*((c.*sigma-qe.*qeprime.*sigmae./kappa).^2)./c; 

r=rr./Ge;% rewritten, checked
rrr=rho+gamma.*i-gamma.*(qe.*qeprime.*(musigmae)+0.5*(qe.*qeprime2+qeprime.^2).*sigmae.^2)./c./kappa...
    -0.5*gamma.*(1+gamma).*(c.*sigma-qe.*qeprime.*sigmae./kappa).^2./c.^2;
SR=maversion*sigma*ttheta.*(we-espan.*weprime)./(we-espan.*m.*ttheta.*weprime); %checked
sigmaRh=(pe*sigma+sigmae.*peprime)./pe; % correct
sigmaRk=(qe*sigma+sigmae.*qeprime)./qe; % correct
rp_h=(peprime.*musigmae+(1/2).*peprime2.*sigmae.^2+((1-phi)./phi).*c)./pe+i-r; % checked
rp_k=(qeprime.*musigmae+(1/2).*qeprime2.*sigmae.^2+A)./qe-delta-r; %checked
muc=i-((qe.*qeprime.*musigmae+(1/2).*(qe.*qeprime2+qeprime.^2).*sigmae.^2)...
    ./(c.*kappa));
sigmac=abs((c.*sigma-qe.*qeprime.*sigmae./kappa)./c);
Doverp=(1-phi)/phi*c./pe;
%SR=muRminusr./sigmaR;
%musigmae2=espan.*(m.^2.*ttheta.^2.*(sigma+sigmae.*weprime./we).^2-eta-i);
qerror=((qeprime.*musigmae+0.5*qeprime2.*sigmae.^2+A)./qe-delta-r)./((sigmae.*qeprime+qe.*sigma)./qe)-SR;
perror=((peprime.*musigmae+0.5*peprime2.*sigmae.^2+(1-phi)/phi.*c)./pe+i-r)./((sigmae.*peprime+pe.*sigma)./pe)-SR;
intequity=min(espan,(1-lambda).*we);
% calculating constraint point and drift
tt=K;ttt=1;
newtt=floor((tt+ttt)/2);
err=we(newtt)-espan(newtt)/(1-lambda);
while tt-ttt>1
    if err>0;
        ttt=newtt;
        newtt=floor((newtt+tt)/2);
    else
        tt=newtt;
        newtt=floor((newtt+ttt)/2);
    end
    err=we(newtt)-espan(newtt)/(1-lambda);    
end
tc=tt; x=(espan(tt)-(1-lambda)*we(tt))/(espan(tt)-espan(ttt)+(1-lambda)*(we(ttt)-we(tt)));
criticalc=espan(tt)+x*(espan(ttt)-espan(tt));
wc=criticalc/(1-lambda);
mue=musigmae-sigma.*sigmae;

%% Output Figures
figure(3)
subplot(531)
plot(espan,pe); title('p(e)'); 
subplot(532); plot(espan,qe); title('q(e)')
subplot(533); plot(espan,SR); title('Sharpe Ratio')
subplot(534); plot(espan,r,espan,rrr); title('Interest rate')
subplot(535); plot(espan,sigmaRk); title('Return Volatility of capital')
subplot(536); plot(espan,rp_k); title('Risk Premium of capital')
subplot(537); plot(espan,sigmaRh); title('Return Volatility of housing')
subplot(538); plot(espan,rp_h); title('Risk Premium of housing')
subplot(539); plot(espan,mueper); title('growth rate of e')
subplot(5,3,10); plot(espan,sigmae./espan); title('percentage volatility of e')
subplot(5,3,11); plot(espan,inv); title('gross investment including adjustment')
subplot(5,3,12); plot(espan,sigmac); title('consumption growth vol')
subplot(5,3,13); plot(espan,c); title('consumption')
subplot(5,3,14); plot(espan,musigmae); title('musigma')
subplot(5,3,15); plot(espan,sigmae./espan); title('e volatility')

end

%% **********  Section: Simulation  ********** 
ec_new=0.7;   count=0;
while( count <= 1 )
    ec = ec_new;
    %load solutiontestsigma4.mat;
    baselineSR=0.319181;
    leverage=(we-intequity)./we;
    TT=4;dt=1/TT;%sim   ulation is quarterly
    he=2;
    num=he*200;  % The original one is 2000. This saves a lot of time, but still quite accurate.
    length1=he*250;length=he*1000;KKlbar=0.01;
    interval=12; %monthly
    %delta=yb/K-0.001;n_b=zeros(K,1);p=zeros(K,1);position=zeros(K,K);Sharpe_position
    %p_bankruptcy=zeros(K,1);p_time=zeros(K,1);deno=zeros(K,1);
    K=500;counter=0;position=zeros(K,1);positionnew=zeros(K,1);deltae=(eubar-elbarnew)/K;
    %J=500;Sharpe_position=zeros(J,1);deltasharpe=gamma/J;
    checkrp=0;
    num_con=0;num_con2=0;n1sigma=0;n2sigma=0;n3sigma=0;n4sigma=0;n5sigma=0;p_time=0;p=0;d=0;averleng=0;
    %we have to adjust the critical value to make sure 33% of distress state.

    %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;std_rf=0;std_rf_con=0;std_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;
    %var_num=4;J=4;
    %e_IRF0=zeros(var_num,J+1,var_num);e_IRF0l=zeros(var_num,J+1,var_num);e_IRF0h=zeros(var_num,J+1,var_num);
    %e_IRF1=zeros(var_num,J+1,var_num);e_IRF1l=zeros(var_num,J+1,var_num);e_IRF1h=zeros(var_num,J+1,var_num);

    cov_equity_equity_con=0;cov_equity_equity_un=0;cov_equity_equity=0;
    cov_equity_inv_con=0;cov_equity_inv_un=0;cov_equity_inv=0;
    cov_equity_sharpe_con=0;cov_equity_sharpe_un=0;cov_equity_sharpe=0;
     cov_equity_consumption_con=0;cov_equity_consumption_un=0;cov_equity_consumption=0;
    cov_equity_landprice_con=0;cov_equity_landprice_un=0;cov_equity_landprice=0;
    cov_sharpe_sharpe_con=0;cov_sharpe_sharpe_un=0;cov_sharpe_sharpe=0;
    cov_sharpe_inv_con=0;cov_sharpe_inv_un=0;cov_sharpe_inv=0;
    cov_sharpe_consumption_con=0;cov_sharpe_consumption_un=0;cov_sharpe_consumption=0;
    cov_sharpe_landprice_con=0;cov_sharpe_landprice_un=0;cov_sharpe_landprice=0;
    cov_inv_inv_con=0;cov_inv_inv_un=0;cov_inv_inv=0;
    cov_inv_landprice_con=0;cov_inv_landprice_un=0;cov_inv_landprice=0;
    cov_inv_consumption_con=0;cov_inv_consumption_un=0;cov_inv_consumption=0;
    cov_consumption_consumption_con=0;cov_consumption_consumption_un=0;cov_consumption_consumption=0;
    cov_consumption_landprice_con=0;cov_consumption_landprice_un=0;cov_consumption_landprice=0;
    cov_landprice_landprice_con=0;cov_landprice_landprice_un=0;cov_landprice_landprice=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;
    mean_equity_growth=0;mean_inv_growth=0;mean_sharpe=0;mean_consumption_growth=0;mean_landprice_growth=0;

    xx=numel(qeprime2);num_con=0;num_un=0;xxx=0;
    tic;
    for x=1:num
        display(  [  'Progress: ' , num2str(100*x/num), '%. Time used: ',  num2str(toc/60,2), ' minutes'   ]   )
        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);
        consum_scaled=zeros(TT*length,1);landprice=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);
        %noise=zeros(TT*length,1);sigmastock=zeros(TT*length,1);
        %portiondele=zeros(TT*length,1);alphaI=zeros(TT*length,1);

        %started coding here
        estart=ec;t=1;KKstart=1;
            while t<TT*length1
            tt=min(floor(1+(estart-elbarnew)/de),xx-1);ttt=tt+1;weight=(estart-(tt-1)*de-elbarnew)/de;
            mu_e=weight*mue(ttt)+(1-weight)*mue(tt);
            sigma_e=weight*sigmae(ttt)+(1-weight)*sigmae(tt);
            dwienerp=sqrt(dt)*randn;
            %inv=weight*i(ttt)+(1-weight)*i(tt);
            %awienerp=awienerp+dwienerp;
            estart=estart+mu_e*dt+sigma_e*dwienerp;

            estart=max(min(estart,eubar),elbarnew);
            t=t+1;

            end
        % starts simulating.    
        ee(1)=estart;indicator=ones(TT*length,1);KK(1)=1;equity(1)=estart;
        t=2;
        %while (hitting==0) & (t<TT*length+2)

            while t<TT*length+2
            %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*mueper(ttt)+(1-weight)*mueper(tt);%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;    
            %if rp(t-1)>=0.05;
            if sharpe(t-1)>=baselineSR*1;
                n1sigma=n1sigma+1;
            end
            %if rp(t-1)>=0.075;
            if sharpe(t-1)>=baselineSR*2; 
            n2sigma=n2sigma+1;
            end
            if sharpe(t-1)>=baselineSR*4;%0.1;
                n3sigma=n3sigma+1;
            end
            if sharpe(t-1)>=baselineSR*8;%rpave*4;%0.2;
                n4sigma=n4sigma+1;
            end
            if sharpe(t-1)>=baselineSR*16;%rpave*4;%0.2;
                n5sigma=n5sigma+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;

           j=max(min(K,floor((ee(t-1))/deltae)+1),1); %NOTE: I record wi along the simulation, not the end of simulation.
        position(j)=position(j)+1;

            t=t+1;
            end

            ee=ee(1:TT*length);
        indi=(ee<=ec);
        indi_growth=indi(1:TT*length-1);
        indi2=(ee<=criticalc);

        if(sum(indi_growth)<=1)   % Newly added by Wenhao to avoid breakdown of the algorithm. As eta becomes small, it is possible to have indi_growht<=1
            continue;
        end
        inv_growth(1:TT*length-1)=log(invest(2:TT*length))-log(invest(1:TT*length-1)); %Investment level
        equity_growth(1:TT*length-1)=log(equity(2:TT*length))-log(equity(1:TT*length-1)); %equity growth
        consumption_growth(1:TT*length-1)=log(consumption(2:TT*length))-log(consumption(1:TT*length-1)); %consumption growth
        landprice_growth(1:TT*length-1)=log(landprice(2:TT*length))-log(landprice(1:TT*length-1)); %landprice growth

        num_con=num_con+sum(indi);num_con2=num_con2+sum(indi2);
        e_ee=e_ee+mean(ee);e_median=e_median+median(ee);e_ee_con=e_ee_con+sum(ee.*indi)/sum(indi);e_ee_un=e_ee_un+sum(ee.*(1-indi))/(TT*length-sum(indi));
        std_ee=std_ee+std(ee);std_ee_con=std_ee_con+std(nonzeros(ee.*indi));std_ee_un=std_ee_un+std(nonzeros(ee.*(1-indi)));

        e_cgrowth=e_cgrowth+mean(cgrowth);e_cgrowth_con=e_cgrowth_con+sum(cgrowth.*indi)/sum(indi);e_cgrowth_un=e_cgrowth_un+sum(cgrowth.*(1-indi))/(TT*length-sum(indi));
        e_cvol=e_cvol+mean(cvol);e_cvol_con=e_cvol_con+sum(cvol.*indi)/sum(indi);e_cvol_un=e_cvol_un+sum(cvol.*(1-indi))/(TT*length-sum(indi));
        e_investment=e_investment+mean(investment);e_investment_con=e_investment_con+sum(investment.*indi)/sum(indi);e_investment_un=e_investment_un+sum(investment.*(1-indi))/(TT*length-sum(indi));
        e_investmentvol=e_investmentvol+mean(investmentvol);e_investmentvol_con=e_investmentvol_con+sum(investmentvol.*indi)/sum(indi);e_investmentvol_un=e_investmentvol_un+sum(investmentvol.*(1-indi))/(TT*length-sum(indi));

        e_pvol=e_pvol+mean(pricevol);e_pvol_con=e_pvol_con+sum(pricevol.*indi)/sum(indi);e_pvol_un=e_pvol_un+sum(pricevol.*(1-indi))/(TT*length-sum(indi));

        %e_invvol=e_invvol+std(invgrowth);e_invvol_con=e_invvol_con+std(nonzeros(invvolgrowth.*indi));e_invvol_un=e_invvol_un+std(nonzeros(invvolgrowth.*(1-indi)))

        e_rentalpriceratio=e_rentalpriceratio+mean(rentalpriceratio);e_rentalpriceratio_con=e_rentalpriceratio_con+sum(rentalpriceratio.*indi)/sum(indi);...
            e_rentalpriceratio_un=e_rentalpriceratio_un+sum(rentalpriceratio.*(1-indi))/(TT*length-sum(indi));

        e_alphak=e_alphak+mean(alphak);e_alphak_con=e_alphak_con+sum(alphak.*indi)/sum(indi);e_alphak_un=e_alphak_un+sum(alphak.*(1-indi))/(TT*length-sum(indi));
        e_alphah=e_alphah+mean(alphah);e_alphah_con=e_alphah_con+sum(alphah.*indi)/sum(indi);e_alphah_un=e_alphah_un+sum(alphah.*(1-indi))/(TT*length-sum(indi));
        e_rpk=e_rpk+mean(rpk);e_rpk_con=e_rpk_con+sum(rpk.*indi)/sum(indi);e_rpk_un=e_rpk_un+sum(rpk.*(1-indi))/(TT*length-sum(indi));
        e_rph=e_rph+mean(rph);e_rph_con=e_rph_con+sum(rph.*indi)/sum(indi);e_rph_un=e_rph_un+sum(rph.*(1-indi))/(TT*length-sum(indi));

        e_rf=e_rf+mean(rf);e_rf_con=e_rf_con+sum(rf.*indi)/sum(indi);e_rf_un=e_rf_un+sum(rf.*(1-indi))/(TT*length-sum(indi));

        e_sigmak=e_sigmak+mean(sigmak);e_sigmak_con=e_sigmak_con+sum(sigmak.*indi)/sum(indi);e_sigmak_un=e_sigmak_un+sum(sigmak.*(1-indi))/(TT*length-sum(indi));
        e_sigmah=e_sigmah+mean(sigmah);e_sigmah_con=e_sigmah_con+sum(sigmah.*indi)/sum(indi);e_sigmah_un=e_sigmah_un+sum(sigmah.*(1-indi))/(TT*length-sum(indi));
        e_sharpe=e_sharpe+mean(sharpe);e_sharpe_con=e_sharpe_con+sum(sharpe.*indi)/sum(indi);e_sharpe_un=e_sharpe_un+sum(sharpe.*(1-indi))/(TT*length-sum(indi));

        std_rf=std_rf+sqrt(TT)*std(rf);std_rf_con=std_rf_con+sqrt(TT)*std(nonzeros(rf.*indi));std_rf_un=std_rf_un+sqrt(TT)*std(nonzeros(rf.*(1-indi)));

        std_serk=std_serk+std(serk);std_serk_un=std_serk_un+std(nonzeros(serk.*(1-indi)));std_serk_con=std_serk_con+std(nonzeros(serk.*indi));
        std_serh=std_serh+std(serh);std_serh_un=std_serh_un+std(nonzeros(serh.*(1-indi)));std_serh_con=std_serh_con+std(nonzeros(serh.*indi));
        std_qq=std_qq+std(qq);std_qq_con=std_qq_con+std(nonzeros(qq.*indi));std_qq_un=std_qq_un+std(nonzeros(qq.*(1-indi)));


        e_theta=e_theta+mean(theta);e_theta_con=e_theta_con+sum(theta.*indi)/sum(indi);e_theta_un=e_theta_un+sum(theta.*(1-indi))/(TT*length-sum(indi));
        e_ppoverqq=e_ppoverqq+mean(ppoverqq);
        e_qq=e_qq+mean(qq);e_qq_con=e_qq_con+sum(qq.*indi)/sum(indi);e_qq_un=e_qq_un+sum(qq.*(1-indi))/(TT*length-sum(indi));

        sharpenew=sharpe(2:TT*length);
        vector_con=[equity_growth(indi_growth==1),inv_growth(indi_growth==1),sharpenew(indi_growth==1),consumption_growth(indi_growth==1),landprice_growth(indi_growth==1)];
        covariance_con=cov(vector_con);mean_con=mean(vector_con);
        cov_equity_equity_con=cov_equity_equity_con+covariance_con(1,1);
        cov_equity_inv_con=cov_equity_inv_con+covariance_con(1,2);
        cov_equity_sharpe_con=cov_equity_sharpe_con+covariance_con(1,3);
        cov_equity_consumption_con=cov_equity_consumption_con+covariance_con(1,4);
        cov_equity_landprice_con=cov_equity_landprice_con+covariance_con(1,5);
        cov_inv_inv_con=cov_inv_inv_con+covariance_con(2,2);
        cov_sharpe_inv_con=cov_sharpe_inv_con+covariance_con(2,3);
        cov_inv_consumption_con=cov_inv_consumption_con+covariance_con(2,4);
        cov_inv_landprice_con=cov_inv_landprice_con+covariance_con(2,5);
        cov_sharpe_sharpe_con=cov_sharpe_sharpe_con+covariance_con(3,3);
        cov_sharpe_consumption_con=cov_sharpe_consumption_con+covariance_con(3,4);
        cov_sharpe_landprice_con=cov_sharpe_landprice_con+covariance_con(3,5);
        cov_consumption_consumption_con=cov_consumption_consumption_con+covariance_con(4,4);
        cov_consumption_landprice_con=cov_consumption_landprice_con+covariance_con(4,5);
        cov_landprice_landprice_con=cov_landprice_landprice_con+covariance_con(5,5);

        mean_equity_growth_con=mean_equity_growth_con+mean_con(1);
        mean_inv_growth_con=mean_inv_growth_con+mean_con(2);
        mean_sharpe_con=mean_sharpe_con+mean_con(3);
        mean_consumption_growth_con=mean_consumption_growth_con+mean_con(4);
        mean_landprice_growth_con=mean_landprice_growth_con+mean_con(5);

        vector_un=[equity_growth(indi_growth==0),inv_growth(indi_growth==0),sharpenew(indi_growth==0),consumption_growth(indi_growth==0),landprice_growth(indi_growth==0)];
        covariance_un=cov(vector_un);mean_un=mean(vector_un);
        cov_equity_equity_un=cov_equity_equity_un+covariance_un(1,1);
        cov_equity_inv_un=cov_equity_inv_un+covariance_un(1,2);
        cov_equity_sharpe_un=cov_equity_sharpe_un+covariance_un(1,3);
        cov_equity_consumption_un=cov_equity_consumption_un+covariance_un(1,4);
        cov_equity_landprice_un=cov_equity_landprice_un+covariance_un(1,5);
        cov_inv_inv_un=cov_inv_inv_un+covariance_un(2,2);
        cov_sharpe_inv_un=cov_sharpe_inv_un+covariance_un(2,3);
        cov_inv_consumption_un=cov_inv_consumption_un+covariance_un(2,4);
        cov_inv_landprice_un=cov_inv_landprice_un+covariance_un(2,5);
        cov_sharpe_sharpe_un=cov_sharpe_sharpe_un+covariance_un(3,3);
        cov_sharpe_consumption_un=cov_sharpe_consumption_un+covariance_un(3,4);
        cov_sharpe_landprice_un=cov_sharpe_landprice_un+covariance_un(3,5);
        cov_consumption_consumption_un=cov_consumption_consumption_un+covariance_un(4,4);
        cov_consumption_landprice_un=cov_consumption_landprice_un+covariance_un(4,5);
        cov_landprice_landprice_un=cov_landprice_landprice_un+covariance_un(5,5);

        mean_equity_growth_un=mean_equity_growth_un+mean_un(1);
        mean_inv_growth_un=mean_inv_growth_un+mean_un(2);
        mean_sharpe_un=mean_sharpe_un+mean_un(3);
        mean_consumption_growth_un=mean_consumption_growth_un+mean_un(4);
        mean_landprice_growth_un=mean_landprice_growth_un+mean_un(5);

        vectortotal=[equity_growth(1:TT*length-1),inv_growth(1:TT*length-1),sharpenew,consumption_growth(1:TT*length-1),landprice_growth(1:TT*length-1)];
        totalcovariance=cov(vectortotal);totalmean=mean(vectortotal);
        cov_equity_equity=cov_equity_equity+totalcovariance(1,1);
        cov_equity_inv=cov_equity_inv+totalcovariance(1,2);
        cov_equity_sharpe=cov_equity_sharpe+totalcovariance(1,3);
        cov_equity_consumption=cov_equity_consumption+totalcovariance(1,4);
        cov_equity_landprice=cov_equity_landprice+totalcovariance(1,5);
        cov_inv_inv=cov_inv_inv+totalcovariance(2,2);
        cov_sharpe_inv=cov_sharpe_inv+totalcovariance(2,3);
        cov_inv_consumption=cov_inv_consumption+totalcovariance(2,4);
        cov_inv_landprice=cov_inv_landprice+totalcovariance(2,5);
        cov_sharpe_sharpe=cov_sharpe_sharpe+totalcovariance(3,3);
        cov_sharpe_consumption=cov_sharpe_consumption+totalcovariance(3,4);
        cov_sharpe_landprice=cov_sharpe_landprice+totalcovariance(3,5);
        cov_consumption_consumption=cov_consumption_consumption+totalcovariance(4,4);
        cov_consumption_landprice=cov_consumption_landprice+totalcovariance(4,5);
        cov_landprice_landprice=cov_landprice_landprice+totalcovariance(5,5);

        mean_equity_growth=mean_equity_growth+totalmean(1);
        mean_inv_growth=mean_inv_growth+totalmean(2);
        mean_sharpe=mean_sharpe+totalmean(3);
        mean_consumption_growth=mean_consumption_growth_un+totalmean(4);
        mean_landprice_growth=mean_landprice_growth+totalmean(5);
    end

    n1sigma=n1sigma/num/TT/length;n2sigma=n2sigma/num/TT/length;
    n3sigma=n3sigma/num/TT/length;n4sigma=n4sigma/num/TT/length;n5sigma=n5sigma/num/TT/length;

    e_ee=e_ee/num;e_median=e_median/num;e_ee_con=e_ee_con/num;e_ee_un=e_ee_un/num;
    std_ee=std_ee/num;std_ee_con=std_ee_con/num;std_ee_un=std_ee_un/num;
    %e_invvol=e_invvol/num;e_invvol_con=e_invvol_con/num;e_invvol_un/num;

    e_cgrowth=e_cgrowth/num;e_cgrowth_con=e_cgrowth_con/num;e_cgrowth_un=e_cgrowth_un/num;
    e_cvol=e_cvol/num;e_cvol_con=e_cvol_con/num;e_cvol_un=e_cvol_un/num;
    e_investment=e_investment/num;e_investment_con=e_investment_con/num;e_investment_un=e_investment_un/num;
    e_investmentvol=e_investmentvol/num;e_investmentvol_con=e_investmentvol_con/num;e_investmentvol_un=e_investmentvol_un/num;
    e_pvol=e_pvol/num;e_pvol_con=e_pvol_con/num;e_pvol_un=e_pvol_un/num;
    e_rentalpriceratio=e_rentalpriceratio/num;e_rentalpriceratio_con=e_rentalpriceratio_con/num;...
        e_rentalpriceratio_un=e_rentalpriceratio_un/num;

    e_alphak=e_alphak/num;e_alphak_con=e_alphak_con/num;e_alphak_un=e_alphak_un/num;
    e_alphah=e_alphah/num;e_alphah_con=e_alphah_con/num;e_alphah_un=e_alphah_un/num;
    e_rpk=e_rpk/num;e_rpk_con=e_rpk_con/num;e_rpk_un=e_rpk_un/num;
    e_rph=e_rph/num;e_rph_con=e_rph_con/num;e_rph_un=e_rph_un/num;

    e_rf=e_rf/num;e_rf_con=e_rf_con/num;e_rf_un=e_rf_un/num;
    e_sigmak=e_sigmak/num;e_sigmak_con=e_sigmak_con/num;e_sigmak_un=e_sigmak_un/num;
    e_sigmah=e_sigmah/num;e_sigmah_con=e_sigmah_con/num;e_sigmah_un=e_sigmah_un/num;
    e_sharpe=e_sharpe/num;e_sharpe_con=e_sharpe_con/num;e_sharpe_un=e_sharpe_un/num;

    std_rf=std_rf/num;std_rf_un=std_rf_un/num;std_rf_con=std_rf_con/num;
    std_serk=std_serk/num/sqrt(dt);std_serk_con=std_serk_con/num/sqrt(dt);std_serk_un=std_serk_un/num/sqrt(dt);
    std_serh=std_serh/num/sqrt(dt);std_serh_con=std_serh_con/num/sqrt(dt);std_serh_un=std_serh_un/num/sqrt(dt);
    std_qq=std_qq/num/sqrt(dt);std_qq_con=std_qq_con/num/sqrt(dt);std_qq_un=std_qq_un/num/sqrt(dt);

    e_theta=e_theta/num;e_theta_con=e_theta_con/num;e_theta_un=e_theta_un/num;
    e_ppoverqq=e_ppoverqq/num;
    e_qq=e_qq/num;e_qq_un=e_qq_un/num;e_qq_con=e_qq_con/num;

    cov_equity_equity_con=cov_equity_equity_con/num;cov_equity_equity_un=cov_equity_equity_un/num;cov_equity_equity=cov_equity_equity/num;
    cov_equity_inv_con=cov_equity_inv_con/num;cov_equity_inv_un=cov_equity_inv_un/num;cov_equity_inv=cov_equity_inv/num;
    cov_equity_sharpe_con=cov_equity_sharpe_con/num;cov_equity_sharpe_un=cov_equity_sharpe_un/num;cov_equity_sharpe=cov_equity_sharpe/num;
    cov_equity_consumption_con=cov_equity_consumption_con/num;cov_equity_consumption_un=cov_equity_consumption_un/num;cov_equity_consumption=cov_equity_consumption/num;
    cov_equity_landprice_con=cov_equity_landprice_con/num;cov_equity_landprice_un=cov_equity_landprice_un/num;cov_equity_landprice=cov_equity_landprice/num;
    cov_sharpe_sharpe_con=cov_sharpe_sharpe_con/num;cov_sharpe_sharpe_un=cov_sharpe_sharpe_un/num;cov_sharpe_sharpe=cov_sharpe_sharpe/num;
    cov_sharpe_inv_con=cov_sharpe_inv_con/num;cov_sharpe_inv_un=cov_sharpe_inv_un/num;cov_sharpe_inv=cov_sharpe_inv/num;
    cov_sharpe_consumption_con=cov_sharpe_consumption_con/num;cov_sharpe_consumption_un=cov_sharpe_consumption_un/num;cov_sharpe_consumption=cov_sharpe_consumption/num;
    cov_sharpe_landprice_con=cov_sharpe_landprice_con/num;cov_sharpe_landprice_un=cov_sharpe_landprice_un/num;cov_sharpe_landprice=cov_sharpe_landprice/num;
    cov_inv_inv_con=cov_inv_inv_con/num;cov_inv_inv_un=cov_inv_inv_un/num;cov_inv_inv=cov_inv_inv/num;
    cov_inv_landprice_con=cov_inv_landprice_con/num;cov_inv_landprice_un=cov_inv_landprice_un/num;cov_inv_landprice=cov_inv_landprice/num;
    cov_inv_consumption_con=cov_inv_consumption_con/num;cov_inv_consumption_un=cov_inv_consumption_un/num;cov_inv_consumption=cov_inv_consumption/num;
    cov_consumption_consumption_con=cov_consumption_consumption_con/num;cov_consumption_consumption_un=cov_consumption_consumption_un/num;cov_consumption_consumption=cov_consumption_consumption/num;
    cov_consumption_landprice_con=cov_consumption_landprice_con/num;cov_consumption_landprice_un=cov_consumption_landprice_un/num;cov_consumption_landprice=cov_consumption_landprice/num;
    cov_landprice_landprice_con=cov_landprice_landprice_con/num;cov_landprice_landprice_un=cov_landprice_landprice_un/num;cov_landprice_landprice=cov_landprice_landprice/num;

    mean_equity_growth_un=mean_equity_growth_un/num;mean_equity_growth_con=mean_equity_growth_con/num;mean_equity_growth=mean_equity_growth/num;
    mean_inv_growth_un=mean_inv_growth_un/num;mean_inv_growth_con=mean_inv_growth_con/num;mean_inv_growth=mean_inv_growth/num;
    mean_sharpe_un=mean_sharpe_un/num;mean_sharpe_con=mean_sharpe_con/num;mean_sharpe=mean_sharpe/num;
    mean_consumption_growth_un=mean_consumption_growth_un/num;mean_consumption_growth_con=mean_consumption_growth_con/num;mean_consumption_growth=mean_consumption_growth/num;
    mean_landprice_growth_un=mean_landprice_growth_un/num;mean_landprice_growth_con=mean_landprice_growth_con/num;mean_landprice_growth=mean_landprice_growth/num;

    num_con=num_con/num/TT/length;num_con2=num_con2/num/TT/length;

    p_over_pq_un =   mean(pp(indi_growth==0)) / mean(  pp(indi_growth==0) + qq(indi_growth==0) ) ;
    avg_lvg_un = mean(ttheta(indi_growth==0)) ;
    prob_distress  = num_con;   % Measured by e below ec
    prob_crisis  = num_con2;  % Measured by the fraction of time spent in the constrained region.
    vol_landprice_growth = sqrt(cov_landprice_landprice*4);
    vol_landprice_growth_con = sqrt(cov_landprice_landprice_con*4);
    mean_output_growth_un = sqrt(  cov_consumption_consumption_un*4 + 2 * cov_inv_consumption_un*4 + cov_inv_inv_un*4  );
    mean_output_growth = sqrt(  cov_consumption_consumption*4 + 2 * cov_inv_consumption*4 + cov_inv_inv*4  );
 
    %main_output = [ elbarnew; qlbarnew;  lambda; A; phi; gamma; etah;  tau ;  avg_lvg_un;  e_investment_un;  p_over_pq_un;   std_rf_un; prob_crisis; vol_landprice_growth; e_rf   ] ;
    main_output = [  etah;   ec;  lambda;  phi; gamma; tau ;  uppersharpe;   ...         % parameters
        prob_crisis;  prob_distress;    avg_lvg_un;  e_investment_un;  p_over_pq_un;   std_rf_un; vol_landprice_growth;     ...  % Key moments
        mean_sharpe_un ; mean_sharpe;  sqrt(cov_consumption_consumption_un*4);   sqrt(cov_inv_inv_un*4) ;    sqrt(cov_equity_equity_un*4);     sqrt(cov_equity_equity_con*4);    ...
        elbarnew; qlbarnew; wlbarnew; mean(ee);  pe(1);   peprime(1);  peprime(1)/pe(1);     err_to_store;  mis_to_store ] ;   % Other output variables

    table4_output = 100*[  sqrt(cov_equity_equity_con*4),  sqrt(cov_inv_inv_con*4),  sqrt(cov_consumption_consumption_con*4),  sqrt(cov_landprice_landprice_con*4),  sqrt(cov_sharpe_sharpe_con*4), ...
        (cov_equity_inv_con*4), (cov_equity_consumption_con*4),  (cov_equity_landprice_con*4),  (cov_equity_sharpe_con*4)     ...
         sqrt(cov_equity_equity_un*4),  sqrt(cov_inv_inv_un*4),  sqrt(cov_consumption_consumption_un*4),  sqrt(cov_landprice_landprice_un*4),  sqrt(cov_sharpe_sharpe_un*4), ...
        (cov_equity_inv_un*4), (cov_equity_consumption_un*4),  (cov_equity_landprice_un*4),  (cov_equity_sharpe_un*4)   ...
        ]' ;

    display(  ['Probability of distress: ', num2str(prob_distress,3) ] )
    display(  ['Probability of crisis: ', num2str(prob_crisis,3) ] )
    display(  ['33% quantile of e distribution: ', num2str(quantile( ee, 0.33 ),3) ] )
    
    ec_new = quantile( ee, 0.33 );
    count = count +1 ;
end
load chirp; sound(y,Fs); 

csvwrite('simulationresult_KeyResults_temp.csv', main_output);   % Key results
csvwrite('table4_result_maversion=1.8.csv', table4_output);   % Key results




