%% set a seed for reproducability of the random numbers (otherwise they will
% be different every time you run the .m file again)
seed = 1234;
rng(seed);
    
% generate 500 obs from a normal distribution with mean = 1, s.d. = 0.01
log_Z = 1 + 0.01.*randn(501,1);

% share of capital alpha and discount factor (beta)
a = 0.35;
b = 0.98;

% initial capital, consumption and output values at time 0
kt(1,1) = 0.9;
ct(1,1) = (1-a*b)*exp(log_Z(1,1))*kt(1,1)^a;
yt(1,1) = exp(log_Z(1,1))*kt(1,1)^a;
time = 0;

% apply the law of motion of capital for each period from 1 to 500
for t = 1:1:500 
       time(t+1,1) = time(t,1) + 1;
       kt(t+1,1) = a*b*exp(log_Z(t,1))*kt(t,1)^a;
       ct(t+1,1) = (1-a*b)*exp(log_Z(t+1,1))*kt(t+1,1)^a;
       yt(t+1,1) = exp(log_Z(t+1,1))*kt(t+1,1)^a;
end

% plot optimal capital, consumption and output paths over time
subplot(3,1,1), plot(time,kt)
xlabel('Time (t)')
ylabel('Capital (K_t)')
subplot(3,1,2), plot(time,ct)
xlabel('Time')
ylabel('Consumption (c_t)')
subplot(3,1,3), plot(time,yt)
xlabel('Time')
ylabel('Output (y_t)')

%% generate logarithm of output and consumption
logy = log(yt);
logc = log(ct);

% display variance of output and consumption
(std(logy))^2
(std(logc))^2

% display standard deviation of output and consumption
std(logy)
std(logc)