%grid.m
% discretizes the process y_t = rho * y_t-1 + sigma_epsilon * epsilon_t,
% where epsilon_t follows a standard normal distribution.
% Output:
% ygrid is the grid of y_t values
% pai transition probability matrix,
clear all
rho = 0.931654648380119;
sigma_epsilon = 0.036960096814455;
W = 4.2; % width of y grid around its mean (0) in terms of standrad devs.
N1 = 200; % initial number of grid points.
% May be reduced if some points are never visited.
T = 1e7; % Length of the time series of simulated values of output
% (used for estimating the transition probability matrix)
% set randomization seed
% randn('state',0);
rng('default')
stdy= sigma_epsilon/sqrt(1-rho^2); % unconditional variance of y_t
UB = W*stdy; % highest value of y grid
ygrid = linspace(-UB,UB,N1); % linearly spaced grid for y
PAI = zeros(N1); % initialize of transition probability matrix
% Simulate time series for log(gdp_traded)
% initialize simulation
[~,i0] = min(abs(ygrid));
y0 = ygrid(i0); % starting point = value of y in middle of grid
fprintf('Creating probability matrix...\n');
% Drawing
for t=2:T
y1 = y0*rho + randn*sigma_epsilon;
[~,i1] = min(abs(y1-ygrid));
% assign draw to transition probability matrix
PAI(i0,i1) = PAI(i0,i1) + 1;
y0 = y1;
i0 = i1;
end
fprintf('Done!\n');
% eliminate all rows and columns with all elements equal to zero
pai = PAI(sum(PAI,2)~=0,sum(PAI,1)~=0);
pai = pai ./ repmat(sum(pai,2),[1,size(pai,2)]);
keep = find(sum(PAI,2)~=0);
ygrid = ygrid(keep);
ygrid = ygrid(:);
N = length(ygrid);
save grid.mat pai ygrid