%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