%simulateEG.m
% produces time series implied by the Eaton-Gersovitz  model.  

clear all, format compact
load loss_spec
rng('default')

eval(['load eg_' outputloss]) %load solution results 

cpai = cumsum(pai,2); %cumulative prob distribution

%Initial conditions
[~,i] = min(abs(y-1)); %average endowment (index)
j = floor(nd/2); %middle of debt grid (index)
b = 0;  %b=0==>bad standing at the end of prevous period;
        %b=1==>good standing at the end of previous period

%Length of simulation and burning period
T = 1e6;
Tburn = 1e5;

Y = zeros(Tburn+T,1);
Ytilde = zeros(Tburn+T,1);
YA = zeros(Tburn+T,1);
B = zeros(Tburn+T,1);
D = zeros(Tburn+T,1);
C = zeros(Tburn+T,1);
Q = zeros(Tburn+T,1);
F = zeros(Tburn+T,1);
TAU = zeros(Tburn+T,1);
STATE = zeros(Tburn+T,1);
R = zeros(Tburn+1,1);
PM = zeros(Tburn+1,1);
rstar_annual = ((1+rstar)^4-1)*100; %world interest rate (annualized)

fprintf('Simulating...\n');
for t=1:T+Tburn

    STATE(t,1) = sub2ind([ny nd ],i,j);

    rr = rand; %random number  determining reentry if applicable
    F(t,1) = f(i,j); %F(t)=1 if default event in period t, 0 otherwise
    D(t,1) = d(j); 
    Y(t,1) = y(i);
    YA(t,1) = ca(i,j);  %output in autarky

    if (b==0) & (F(t)==0) ; %choose to continue 
        B(t,1) = 0; %B=0==>good standing at the beginning of current period, after
                    %receving reentry signal if applicable, but before making a  
                    %default decision. . 
                    %B==1==>bad  standing at the beginning of current period, after
                    %receving reentry signal 
        Ytilde(t,1) = y(i);
        C(t,1) = cc(i,j);
        Q(t,1) = q(i,dpix(i,j));
        TAU(t,1) = tauc(i,j)*100;
        jp  = dpix(i,j); %update debt state
    end 

    if (b==0) & (F(t) ==1); %choose to default
        B(t,1) = 0; 
        Ytilde(t,1) = YA(t);
        C(t,1) = YA(t);
        Q(t,1) = qa(i,j);
        TAU(t,1) = 0;
        jp = nd0; 
    end 

    if (b==1) & (rr>theta);  %==> autarky (bad standing and did not get to re-enter)
        B(t,1) = 1; 
        Ytilde(t,1) = YA(t);
        C(t,1) = YA(t);
        Q(t,1) = qa(i,j);
        TAU(t,1) = 0;
        jp = nd0; 
    end 

    if (b==1) & (rr<=theta) ; %reentry after having been in bad standing 
        B(t,1) = 0; 
        Ytilde(t,1) = Y(t);
        C(t,1) = cc(i,nd0);
        Q(t,1) = q(i,dpix(i,nd0));
        TAU(t,1) = tauc(i,nd0)*100;
        jp = dpix(i,nd0);
    end

    R(t,1) = ((1/Q(t,1))^4-1)*100; %country interest rate 
    PM(t,1) = R(t,1)-rstar_annual;%country risk premium

    %update output state
    find(cpai(i,:)>rand); 
    i = ans(1);

    b = B(t) + F(t); 

    %update debt state
    j=jp; 

end %for t=1:T+Tburn
fprintf('Done!\n');

%eliminate burning period
STATE = STATE(Tburn+1:end);
Y = Y(Tburn+1:end);
YA = YA(Tburn+1:end);
Ytilde = Ytilde(Tburn+1:end);
B = B(Tburn+1:end);
D = D(Tburn+1:end);
C = C(Tburn+1:end);
Q = Q(Tburn+1:end);
F = F(Tburn+1:end);
TAU = TAU(Tburn+1:end);
PM = PM(Tburn+1:end);
R = R(Tburn+1:end);


Ddistr = NaN(nd,1); %debt distribution
for i=1:nd
    Ddistr(i,1) = mean(D==d(i))*100;
end

eval(['save simu_eg_' outputloss])