clear;
format('compact')

global solF etah kappa A delta rho targetR currente currentK sigma m phi KKK beta lambda criticalc gamma maversion elbarnew de  dt

load simulation_Baseline.mat;
%load simulation_maversion=2.4.mat;   

% Note: ec should be the same as the loaded file
e0=ec; 
dt = 1/4;  % The smallest unit of time is quarter.

% 2008-2010 equitiy changes
Requity=[-0.130311833
-0.144248382
-0.152228461
0.060648304
-0.308922417
-0.203930345
0.36811311
0.292671427
-0.016253934];

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);
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,1);dZ(1)=-0;
for t=1:T
    currente=e(t);currentK=K(t);targetR=R(t);
    dZ(t+1)=fsolve(@myfunequity,dZ(t));
    %dZ(t)=fsolve(@myfunlandprice,0);
    xx=musigma(currente);
    e(t+1)=max(e(t)+xx(1)*dt+xx(2)*dZ(t+1),elbarnew);
    K(t+1)=KKK;
    outcome=ils(e(t+1));
    investment(t+1)=outcome(1)*K(t+1);
    landprice(t+1)=outcome(2)*K(t+1);
    Sharperatio(t+1)=outcome(3);
    equity(t+1)=min((1-lambda)*(outcome(2)+outcome(4)),e(t+1))*K(t+1);
end

equitynorm=equity/equity(1);
landpricenorm=landprice/landprice(1);
investmentnorm=investment/investment(1);
percentageChangesinCapital = sigma*dZ;
display(  ['Cumulative decline in productive capital: ', num2str( 100 - 100*prod( percentageChangesinCapital+1 ) ), '%']  )

output=[equitynorm,investmentnorm,landpricenorm,Sharperatio, percentageChangesinCapital];
xlswrite('Figure3_Simulation.xlsx', output);



%% Added on April 28, 2017. Hidden Leverage
global lambdaprime  
lambdaprime=lambda+0.022;

% With hidden leverage
e_hidden=zeros(T,1);K_hidden=zeros(T,1);investment_hidden=zeros(T,1);landprice_hidden=zeros(T,1);Sharperatio_hidden=zeros(T,1);equity_hidden=zeros(T,1);
e_hidden(1)=e0;  K_hidden(1)=1;

outcome_hidden=ils_hiddenLeverage(e_hidden(1));
investment_hidden(1)=outcome_hidden(1)*K_hidden(1);
landprice_hidden(1)=outcome_hidden(2)*K_hidden(1);
Sharperatio_hidden(1)=outcome_hidden(3);
equity_hidden(1)=min((1-lambdaprime)*(outcome_hidden(2)+outcome_hidden(4)),e_hidden(1))*K_hidden(1);
for t=1:T
    currente=e_hidden(t);currentK=K_hidden(t);
    % Solve the evolution of e
    xx  = musigma_higherleverage(currente);
    mu_e = xx(1);  
    sigma_e = xx(2);
    
    % Solve next period e
    e_hidden(t+1)=max(e_hidden(t)+mu_e*dt+sigma_e*dZ(t+1),elbarnew);
    
    % Solve next period capital K
    KKK = myfunequity_higherleverage(dZ(t+1));
    K_hidden(t+1)=KKK;
    
    % Update other variables
    outcome_hidden=ils_hiddenLeverage(e_hidden(t+1));
    investment_hidden(t+1)=outcome_hidden(1)*K_hidden(t+1);
    landprice_hidden(t+1)=outcome_hidden(2)*K_hidden(t+1);
    Sharperatio_hidden(t+1)=outcome_hidden(3);
    equity_hidden(t+1)=min((1-lambdaprime)*(outcome_hidden(2)+outcome_hidden(4)),e_hidden(t+1))*K_hidden(t+1);
end
equitynorm_hidden=equity_hidden/equity_hidden(1);


% Compare different functions with high and low leverage
riskpremium=qe./we.*rp_k+pe./we.*rp_h;
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);
mueprime=musigmaeprime-sigma.*sigmaeprime;

ttheta=we./espan.*(espan<criticalc)+1./(1-lambda).*(espan>=criticalc);
sigmae=(espan.*we.*sigma.*(m.*ttheta-1))./(we-espan.*m.*ttheta.*weprime); %checked 
musigmae=espan.*(m.*r+m.*ttheta.*(riskpremium-r)-etah-i);
mue=musigmae-sigma.*sigmae;

plot( sigmaeprime - sigmae )
plot( mueprime - mue )
plot( tthetaprime - ttheta )


%% Under higher leverage
load simulation_lambda=0.772.mat;

e_high=zeros(T,1);K_high=zeros(T,1);investment_high=zeros(T,1);landprice_high=zeros(T,1);Sharperatio_high=zeros(T,1);equity_high=zeros(T,1);
e_high(1)=e0;  K_high(1)=1;

outcome_high=ils(e_high(1));
investment_high(1)=outcome_high(1)*K_high(1);
landprice_high(1)=outcome_high(2)*K_high(1);
Sharperatio_high(1)=outcome_high(3);
equity_high(1)=min((1-lambda)*(outcome_high(2)+outcome_high(4)),e_high(1))*K_high(1);
for t=1:T
    currente=e_high(t);currentK=K_high(t);
    % Solve the evolution of e
    xx  = musigma(currente);
    mu_e = xx(1);  
    sigma_e = xx(2);
    
    % Solve next period e
    e_high(t+1)=max(e_high(t)+mu_e*dt+sigma_e*dZ(t+1),elbarnew);
    
    % Solve next period capital K
    myfunequity(dZ(t+1));
    K_high(t+1)=KKK;
    
    % Update other variables
    outcome_high=ils(e_high(t+1));
    investment_high(t+1)=outcome_high(1)*K_high(t+1);
    landprice_high(t+1)=outcome_high(2)*K_high(t+1);
    Sharperatio_high(t+1)=outcome_high(3);
    equity_high(t+1)=min((1-lambda)*(outcome_high(2)+outcome_high(4)),e_high(t+1))*K_high(t+1);
end
equitynorm_high=equity_high/equity_high(1);


Table = [ equitynorm, Sharperatio, equitynorm_hidden, Sharperatio_hidden, equitynorm_high,  Sharperatio_high   ];

xlswrite('ThreeCases_Leverage.xlsx', Table);




