function [Z, P] = tauchen(rho, sigma_eps, N, m)
% TAUCHEN Discretizes an AR(1) process using Tauchen's (1986) method.
%
%   [Z, P] = tauchen(rho, sigma_eps, N, m)
%
%   Inputs:
%       rho       - Persistence parameter of the AR(1) process (0 <= rho < 1)
%       sigma_eps - Standard deviation of the innovation term (epsilon)
%       N         - Number of states in the discretization
%       m         - Number of standard deviations to cover (typically m=3)
%
%   Outputs:
%       Z         - (N x 1) vector of state values (the grid for z)
%       P         - (N x N) transition probability matrix. P(i,j) is prob
%                   of moving from state i to state j.
%
%   The AR(1) process is: z_{t+1} = rho * z_t + epsilon_t
%   where epsilon_t ~ N(0, sigma_eps^2)

    % 1. Calculate the unconditional standard deviation of z
    sigma_z = sigma_eps / sqrt(1 - rho^2);

    % 2. Create the grid Z
    % The grid is centered at 0 and covers +/- m standard deviations
    z_min = -m * sigma_z;
    z_max =  m * sigma_z;
    Z = linspace(z_min, z_max, N)';
    
    % Step size between grid points
    d = Z(2) - Z(1);

    % 3. Compute the transition probability matrix P
    P = zeros(N, N);
    
    % Loop over current state i and future state j
    for i = 1:N
        for j = 1:N
            % Define the boundaries for the interval corresponding to state j
            % Upper bound: Z(j) + d/2
            % Lower bound: Z(j) - d/2
            
            % Adjust boundaries for the first and last states to cover tails
            if j == 1
                % Probability of being in the left tail
                prob = normcdf((Z(j) + d/2 - rho*Z(i)) / sigma_eps);
            elseif j == N
                % Probability of being in the right tail
                prob = 1 - normcdf((Z(j) - d/2 - rho*Z(i)) / sigma_eps);
            else
                % Probability of being in the interval centered at Z(j)
                prob = normcdf((Z(j) + d/2 - rho*Z(i)) / sigma_eps) - ...
                       normcdf((Z(j) - d/2 - rho*Z(i)) / sigma_eps);
            end
            
            P(i, j) = prob;
        end
    end
    
    % Ensure rows sum to 1 (handle minor numerical errors)
    P = P ./ sum(P, 2);

end

