% Code adapted from the work by Mattias Almgren (IIES)
clear all; 
clc;

%% set parameters and starting values
par.rho                 = 0.9;      % AR(1) parameter for income process
par.sigma_varepsilon    = 0.03;     % variance of output shock
ybar                    = 1;        % mean output
par.r                   = 0.1;      % exogeneous interest rate
d_minus1                = ybar/2;   % starting value savings
y_minus1                = ybar;     % starting value income

%% Simulation
% simulate output process once
par.T           = 100;      % number periods to simulate.
par.Tdiscard    = 50;       % number periods to discard.
par.Niter       = 1000;     % number of replications

moments         = NaN(par.Niter,3); % pre-allocate matrix for storage of moments
rng(1234)                           % seed for random number generator. For reproducibility. 
for iter = 1:1000
    % preallocate vectors to store realized outcome
    ytemp       = NaN(par.T,1);     
    ctemp       = NaN(par.T,1);
    dtemp       = NaN(par.T,1);
    tbtemp      = NaN(par.T,1);
    ytemp(1)    = y_minus1;         % manually insert starting value for output
    dtemp(1)    = d_minus1;         % manually insert starting value for savings
    
    for tt = 2:par.T
        ytemp(tt)   = par.rho * ( ytemp(tt-1) - ybar ) + ybar + random('Normal',0,par.sigma_varepsilon,1,1);
        ctemp(tt)   = ybar + (par.r / (1 + par.r - par.rho) ) * ( ytemp(tt) - ybar ) - par.r*dtemp(tt-1);
        dtemp(tt)   = dtemp(tt-1) - ( (1-par.rho) / ( 1 + par.r - par.rho ) ) * ( ytemp(tt) - ybar );
        tbtemp(tt)  = ( (1-par.rho) / (1+par.r-par.rho) ) * (ytemp(tt) - ybar) + par.r * dtemp(tt-1);
    end

    % compute trade-balance-to-output
    tby     = tbtemp./ytemp;

    % calculate growth rates, and discard par.Tdiscard first periods
    gy = log(ytemp) - lagmatrix(log(ytemp),1);
    gy = gy(par.Tdiscard+1:end);

    gc = log(ctemp) - lagmatrix(log(ctemp),1);
    gc = gc(par.Tdiscard+1:end);

    tby = tby(par.Tdiscard+1:end);

    % compute standard deviations of respective growth series
    sigma_gy    = std(gy);
    sigma_gc    = std(gc);
    rho_gytby   = corr(gy,tby);

    % store moments in result matrix
    moments(iter,1) = sigma_gy;
    moments(iter,2) = sigma_gc;
    moments(iter,3) = rho_gytby;

end

average_moments = mean(moments);

disp(['sigma_gy = ',num2str(average_moments(1)),', sigma_gc = ',num2str(average_moments(2)),', rho_gytby = ',num2str(average_moments(3))])