% This code simulates various stochastic processes.
%
% Before running the program, you should set the path of this source file.
% The code starts running after you type the command >> simvar.
%
% After each figure generated, return to the command window
% to inspect the next process to be generated.
%
% This code calls >> normrnd function of the Statistics Toolbox of MATLAB.
%
% (c) Mustafa A. Attar . November 25, 2007

clc
close all
clear all

disp(' ')
disp('y[t+1] = RHO*y[t] + e[t+1]')
disp(' ')
disp('y[0] = 0.00              Initial Value')
disp('e[.] ~ N[0,0.01]         Shock Process')
disp('RHO  = 0.85              Persistence  ')
disp(' ')
disp('> > > > > Press any key to display the simulated AR(1) process')

pause

T = 200;
e = zeros(200,1);
y = zeros(200,1);
rho = .85;
e(2:200,1) = normrnd(0,0.01,199,1);
for t=2:T-1;
    y(t+1,1) = rho*y(t,1) + e(t+1,1);
end;

figure(1)
plot(y(3:200,1))
xlabel('time')
ylabel('A Stationary AR(1) Process')

disp(' ')
disp(' ')
disp('> > > > > Press any key to display the process for RHO = 0.95')
disp(' ')
disp(' ')

pause

T = 200;
e = zeros(200,1);
y = zeros(200,1);
rho = .95;
e(2:200,1) = normrnd(0,0.01,199,1);
for t=2:T-1;
    y(t+1,1) = rho*y(t,1) + e(t+1,1);
end;

figure(2)
plot(y(3:200,1))
xlabel('time')
ylabel('More Persistent AR(1)')

disp(' ')
disp(' ')
disp('> > > > > Press any key to display the process for RHO = 1.00')
disp('          This process is called a "Random Walk"')
disp(' ')
disp(' ')

pause

T = 200;
e = zeros(200,1);
y = zeros(200,1);
rho = 1.00;
e(2:200,1) = normrnd(0,0.01,199,1);
for t=2:T-1;
    y(t+1,1) = rho*y(t,1) + e(t+1,1);
end;

figure(3)
plot(y(3:200,1))
xlabel('time')
ylabel('A Random Walk Process')

disp(' ')
disp(' ')
disp('Random Walk process is difference-stationary. That is, given RHO = 1.00,')
disp('we have y[t+1] - y[t] = e[t+1] which is stationary.')
disp(' ')
disp('The following example simulates a process which is purely trend-stationary.')
disp(' ')
disp('y[t+1] = PHI_0 + PHI_1*(t+1) + PHI_2*(t+1)^2 + e[t+1]')
disp(' ')
disp('In above specification of y[t], the deterministic components are')
disp('   - constant term                  PHI_0')
disp('   - linear trend term        PHI_1*(t+1)')
disp('   - quadratic trend term   PHI_2*(t+1)^2')
disp(' ')
disp(' ')
disp('> > > > > Press any key to display the trend-stationary process.')
disp('          for PHI_0 =  1.00000')
disp('              PHI_1 =  0.00500')
disp('              PHI_2 =  0.00004')
disp('              y[0]  =  0.00000')
disp('              e[.]  ~  N[0,0.1]')

pause

T = 200;
e = zeros(200,1);
y = zeros(200,1);
phi_0 = 1.00;
phi_1 = 0.0050;
phi_2 = 0.00004;
e(2:200,1) = normrnd(0,0.1,199,1);
for t=2:T-1;
    y(t+1,1) = phi_0 + phi_1*(t+1) + phi_2*(t+1)^2 + e(t+1,1);
end;

figure(4)
plot(y(3:200,1))
xlabel('time')
ylabel('A Trend-Stationary Process')

disp(' ')
disp('> > > > > Press any key to continue with vector processes.')
disp(' ')

pause

clc
close all
clear all

disp(' ')
disp(' ')
disp('We can analyze vector processes in a similar way. Consider first the')
disp('following bivariate process:')
disp('   y1[t+1] = f_y1(y1[t],y2[t],e_y1[t+1])')
disp('   y2[t+1] = f_y2(y1[t],y2[t],e_y2[t+1])')
disp('where f_y1 and f_y2 are known functions. In what follows, we assume both')
disp('of these functions are linear. In matrix form, such a bivariate system')
disp('can be written as')
disp('   A*y[t+1] = B*y[t] + C*e[t+1]')
disp('Now, we can specify the model by restricting matrices A, B and C.')
disp('For simplicity, assume A and C are identity matrices of size (2x2)')
disp('and B(i,j) for i,j = 1,2 are known. Then, we have')
disp('   y1[t+1] = B(1,1)*y1[t] + B(1,2)*y2[t] + e_y1[t+1]')
disp('   y2[t+1] = B(2,1)*y1[t] + B(2,2)*y2[t] + e_y2[t+1]')
disp(' ')
disp('For the simulation, assume B(1,1) =  0.9')
disp('                           B(1,2) = -0.2')
disp('                           B(2,1) =  0.1')
disp('                           B(2,2) =  0.7')
disp(' ')
disp('                           e_y1   ~  N[0,1]')
disp('                           e_y2   ~  N[0,1]')
disp(' ')
disp('                           y1[0]  =  0.0')
disp('                           y2[0]  =  0.0')
disp(' ')
disp('> > > > > Press any key now to display the bivariate stationary process.')

pause

T = 200;
E = normrnd(0,1,2,200);
Y = zeros(2,200);
B = zeros(2);
B(1,1) = 0.9; B(1,2) = -0.2; B(2,1) = 0.1; B(2,2) = 0.7;
A = eye(2);
C = A;
E(:,1) = 0;
for t=2:T-1;
    Y(:,t+1) = A^(-1)*B*Y(:,t) + A^(-1)*C*E(:,t+1);
end;

Y = Y';

figure(1)
subplot(2,1,1)
plot(Y(3:200,1))
xlabel('time')
ylabel('y1')
subplot(2,1,2)
plot(Y(3:200,2))
xlabel('time')
ylabel('y2')

disp(' ')
disp(' ')
disp('Note that one can impose contemporaneous relationship between')
disp('y1 and y2 by using A matrix and between e_y1 and e_y2 using')
disp('C matrix. These exercises are left to the reader.')
disp(' ')
disp('We will now discuss a very simple dynamic, stochastic,')
disp('general equilibrium model and simulate it using its state-space')
disp('representation.')
disp(' ')
disp('> > > > > Press any key to continue.')

pause

clc
close all
clear all

disp('We specify a State-Space model as follows:')
disp('       y[t+1] = A*x[t+1]')
disp('       x[t+1] = B*x[t] + C*e[t+1]')
disp('where  y: the vector of control variables')
disp('       x: the vector of state variables')
disp('       e: the vector of shocks (or driving forces)')
disp(' ')
disp('The DSGE model that we will use here is stochastic, one sector, Ramsey')
disp('growth model.')
disp('PREFERENCES')
disp('E_0[log(C_0) + BETTA*log(C_1) + (BETTA^2)*log(C_2) + ...]   (i)')
disp('PRODUCTION')
disp('Y_t = A_t*(K_t)^ALFA                                       (ii)')
disp('CAPITAL ACCUMULATION')
disp('K_t+1 = (1-DELTA)*K_t + I_t                               (iii)')
disp('RESOURCE CONSTRAINT')
disp('Y_t = C_t + I_t                                            (iv)')
disp('EVOLUTION OF TECHNOLOGY')
disp('log(A_t+1) = RHO*log(A_t) + e_t+1)                          (v)')
disp(' ')
disp('Control Vector: y = [Y C I]')
disp('State Vector:   x = [K A]')
disp(' ')
disp('> > > > > Press any key to continue.')

pause

disp('The solution of the Ramsey problem is characterized by the Euler')
disp('equation:')
disp('      C_t+1 / C_t = BETTA*[ALFA*A_t+1*(K_t+1)^(ALFA-1) + (1-DELTA)]')
disp('Euler Equation and Equations (ii)-(v) define a nonlinear system.')
disp('The log-linearization of this system around its non-stochastic')
disp('steady-state transforms the system of y and x from <levels> into')
disp('<percentage deviations from steady-state>.')
disp('Assigning values to the parameters ALFA, BETTA, GAMA and RHO, the')
disp('linear system can be solved to obtain A and B matrices.')
disp('SEE: Schmitt-Grohe and Uribe (2004)')
disp('     in Journal of Economic Dynamics and Control')
disp('> > > > > Press any key to continue.')

pause

disp('PARAMETER VALUES:')
disp('      ALFA  = 0.30')
disp('      BETTA = 0.96')
disp('      DELTA = 0.10')
disp('      RHO   = 0.70')
disp(' ')
disp(' A =  0.3000    1.0000')
disp('      0.4426    0.6743')
disp('     -1.9678    6.1820')
disp(' ')
disp(' B =  0.7032    0.6182')
disp('      0.0000    0.7000')
disp(' ')
disp(' C =  0.0000')
disp('      1.0000')
disp(' ')
disp('> > > > > Press any key to display the simulated series of')
disp('          Output, Consumption and Investment under the assumptions')
disp('          that  1. the shock to technology, e_t, follows N[0,1].')
disp('                2. initial values for all variables are 0.')

pause

T = 100;
e = normrnd(0,1,1,T);
y = zeros(3,T);
x = zeros(2,T);
A = [0.3000 1.0000; 0.4426 0.6743; -1.9678 6.1820];
B = [0.7032 0.6182; 0.0000 0.7000];
C = [0;1];
for t=2:T-1
    x(:,t)=B*x(:,t-1)+C*e(:,t);
    y(:,t)=A*x(:,t);
end
Y = y';
X = x';

figure(1)
subplot(3,1,1)
plot(Y(2:100,1))
title('Deviations from Steady-State Levels')
ylabel('Output')
subplot(3,1,2)
plot(Y(2:100,2))
ylabel('Consumption')
subplot(3,1,3)
plot(Y(2:100,3))
ylabel('Investment')

disp('> > > > > Press any key to display the simulated series of')
disp('          Capital Stock and Technology.')

pause

figure(2)
subplot(2,1,1)
plot(X(3:100,1))
title('Deviations from Steady-State Levels')
ylabel('Capital')
subplot(2,1,2)
plot(X(3:100,2))
ylabel('Technology')

disp(' ')
disp('> > > > > This is the end. Press any key to clear all.')

pause

clc
clear all
close all