Discrete Dynamic Programming

Author

Pablo Winant

Markov Chains

A worker’s employment dynamics obey the stochastic matrix

\[P = \begin{bmatrix} 1-\alpha & \alpha \\ \beta & 1-\beta \end{bmatrix}\]

\[P = \begin{bmatrix} 1-\alpha & ... \\ \beta & ... \end{bmatrix}\]

with \(\alpha\in(0,1)\) and \(\beta\in (0,1)\). First line corresponds to employment, second line to unemployment.

Which is the stationary equilibrium? (choose any value for \(\alpha\) and \(\beta\))

α = 0.123
β = 0.345
P = [1-α α ; β 1-β]
2×2 Matrix{Float64}:
 0.877  0.123
 0.345  0.655
α = 0.123
β = 0.345
γ = 0.678
P =/2 1-α  α/2; 1-β β/2 β/2; γ*0.8 γ*0.2 1-γ]
3×3 Matrix{Float64}:
 0.0615  0.877   0.0615
 0.655   0.1725  0.1725
 0.5424  0.1356  0.322
α = 0.123
β = 0.345
γ = 0.678
P =/2 1-α  α/2; 1-β β/2 β/2; γ*0.8 γ*0.2 1-γ]
P
3×3 Matrix{Float64}:
 0.0615  0.877   0.0615
 0.655   0.1725  0.1725
 0.5424  0.1356  0.322
Pp = P'
3×3 adjoint(::Matrix{Float64}) with eltype Float64:
 0.0615  0.655   0.5424
 0.877   0.1725  0.1356
 0.0615  0.1725  0.322
3-element Vector{Float64}:
 0.41963333333333336
 0.39503333333333335
 0.18533333333333332
μbar = μ0'*P^50 # stupid too many matrix multiplications 
1×3 adjoint(::Vector{Float64}) with eltype Float64:
 0.400406  0.44903  0.150564
μ0 = ones(3)/3 # initial probability vector

for t=1:40
    μ1 = Pp * μ0
    μ0 = μ1
end

μ0
3-element Vector{Float64}:
 0.4004056827203767
 0.4490299796026508
 0.1505643376769728
size(P, 1)
3
# refactor as a function to solve any stochastic matrix by simulation

function solve_by_simulation(P; T=1000, τ_η=1e-10)

    n = size(P, 1)

    Pp = P'

    μ0 = ones(n)/n # initial probability vector

    for t=1:T
        μ1 = Pp * μ0

        # check whether we have a fixed point
        # ϵ = maximum(abs, μ1 - Pp*μ1)
        # not useful because Pp*μ1 will be computed in the next iteration

        # check successsive approximation errors
        η = maximum(abs,μ1 - μ0)

        if η<τ_η
            return μ1
        end

        μ0 = μ1
    end

    error("No convegence")

end
solve_by_simulation (generic function with 1 method)
solve_by_simulation(P)
3-element Vector{Float64}:
 0.40040568284433287
 0.4490299794605776
 0.15056433769508987
### Solve with linear algebra

# let's solve M μ = μ by defining the appropriate M
using LinearAlgebra: I
P' - I
3×3 Matrix{Float64}:
 -0.9385   0.655    0.5424
  0.877   -0.8275   0.1356
  0.0615   0.1725  -0.678
# using least square solver
F = cat( P' - I, ones(1, 3); dims=1)
R = zeros(4)
R[end] = 1.0
F \ R # here it solves a nonsquare system
3-element Vector{Float64}:
 0.4004056828214058
 0.44902997948685536
 0.15056433769173883
# using linear solver
M = cat( P' - I,; dims=1)
# replace last line with ones
M[end,:] .= 1.0
D = zeros(3)
D[end] = 1.0
M \ D
3-element Vector{Float64}:
 0.40040568282140593
 0.44902997948685536
 0.15056433769173877

In the long run, what will the the fraction \(p\) of time spent unemployed? (Denote by \(X_m\) the fraction of dates were one is unemployed)

using Plots


function forecast(P; T=100, μ0=ones(size(P,1))/size(P,1) )

    sim =0]

    Pp = P'

     # initial probability vector

    for t=1:T
        μ1 = Pp * μ0
        push!(sim, μ1)

        μ0 = μ1
    end

    return sim
end
forecast (generic function with 1 method)
forecast(P)
101-element Vector{Vector{Float64}}:
 [0.3333333333333333, 0.3333333333333333, 0.3333333333333333]
 [0.41963333333333336, 0.39503333333333335, 0.18533333333333332]
 [0.3850790833333334, 0.4612928833333334, 0.15362803333333333]
 [0.4091570474883334, 0.4381193397783334, 0.15272361273333335]
 [0.3949686135219009, 0.4551156386456709, 0.14991574783242834]
 [0.4037056146688205, 0.4452234971311626, 0.15107088820001704]
 [0.39839013568273324, 0.45133608975960343, 0.15027377455766347]
 [0.40163462745710504, 0.4476207483073078, 0.1507446242355873]
 [0.3996560039152811, 0.44988911840923734, 0.15045487767548166]
 [0.4008629424500215, 0.4485058697720903, 0.1506311877778883]
 ⋮
 [0.4004056828214061, 0.4490299794868556, 0.15056433769173902]
 [0.40040568282140615, 0.4490299794868555, 0.15056433769173902]
 [0.4004056828214061, 0.4490299794868556, 0.15056433769173902]
 [0.40040568282140615, 0.4490299794868555, 0.15056433769173902]
 [0.4004056828214061, 0.4490299794868556, 0.15056433769173902]
 [0.40040568282140615, 0.4490299794868555, 0.15056433769173902]
 [0.4004056828214061, 0.4490299794868556, 0.15056433769173902]
 [0.40040568282140615, 0.4490299794868555, 0.15056433769173902]
 [0.4004056828214061, 0.4490299794868556, 0.15056433769173902]

Illustrate this convergence by generating a simulated series of length 10000 starting at \(X_0=1\). Plot \(X_m-p\) against \(m\). (Take \(\alpha=\beta=0.1\)).

Basic Asset Pricing model

A financial asset yields dividend \((x_t)\), which follows an AR1. It is evaluated using the stochastic discount factor: \(\rho_{0,t} = \beta^t \exp(y_t)\) where \(\beta<1\) and \(y_t\) is an \(AR1\). The price of the asset is given by \(p_0 = \sum_{t\geq 0} \rho_{0,t} U(x_t)\) where \(U(u)=\exp(u)^{0.5}/{0.5}\). Our goal is to find the pricing function \(p(x,y)\), which yields the price of the asset in any state.

Write down the recursive equation which must be satisfied by \(p\).

\[p_t = U(x_t) + \beta E_t \left[ \frac{e^{y_{t+1}}}{e^{y_t}} p_{t+1} \right]\]

Compute the ergodic distribution of \(x\) and \(y\).

Discretize processes \((x_t)\) and \((y_t)\) using 2 states each. How would you represent the unknown \(p()\)?

Solve for \(p()\) using successive approximations

Solve for \(p()\) by solving a linear system (homework)

Asset replacement (from Compecon)

At the beginning of each year, a manufacturer must decide whether to continue to operate an aging physical asset or replace it with a new one.

An asset that is \(a\) years old yields a profit contribution \(p(a)\) up to \(n\) years, at which point, the asset becomes unsafe and must be replaced by law.

The cost of a new asset is \(c\). What replacement policy maximizes profits?

Calibration: profit \(p(a)=50-2.5a-2.5a^2\). Maximum asset age: 5 years. Asset replacement cost: 75, annual discount factor \(\delta=0.9\).

Define kind of problem, the state space, the actions, the reward function, and the Bellman updating equation

Solve the problem using Value Function Iteration

Solve the problem using Policy Iteration. Compare with VFI.