-
-
Save korenmiklos/3c1b8cdc19fcf27b4da80631882c59e6 to your computer and use it in GitHub Desktop.
| using SpectralKit | |
| using LinearAlgebra | |
| using Optim | |
| using Plots | |
| const maxiter = 10 | |
| mutable struct SaddlePath | |
| basis::SpectralKit.FunctionBasis | |
| coefs::Vector{Float64} | |
| end | |
| costate(sp::SaddlePath) = linear_combination(sp.basis, sp.coefs) | |
| 𝑑costate(sp::SaddlePath) = x -> linear_combination(sp.basis, sp.coefs, 𝑑(x))[1] | |
| struct RamseyParams | |
| α::Float64 # Capital share | |
| β::Float64 # Discount factor | |
| δ::Float64 # Depreciation rate | |
| end | |
| f(k, α) = k^α | |
| 𝑑f(k, α) = α * k^(α-1) | |
| # Steady state | |
| k_ss(p::RamseyParams) = ((1/p.β - 1 + p.δ) / p.α)^(1/(p.α - 1)) | |
| c_ss(p::RamseyParams) = k_ss(p)^p.α - p.δ * k_ss(p) | |
| function saddle_path(p::RamseyParams, window::Float64=0.25, n::Int=5) | |
| k_min = (1 - window) * k_ss(p) | |
| k_max = (1 + window) * k_ss(p) | |
| basis = Chebyshev(EndpointGrid(), n) ∘ BoundedLinear(k_min, k_max) | |
| SaddlePath(basis, zeros(dimension(basis))) | |
| end | |
| rp = RamseyParams(0.36, 0.96, 0.05) | |
| sp = saddle_path(rp, 0.5, 13) | |
| kss = k_ss(rp) | |
| css = c_ss(rp) | |
| cdot(k) = 𝑑f(k, rp.α) - rp.δ + 1 - rp.β | |
| function distance(sp::SaddlePath, f1::Function, f2::Function) | |
| RSS = 0.0 | |
| for point in SpectralKit.grid(sp.basis) | |
| RSS += (f1(point) - f2(point))^2 | |
| end | |
| return RSS / length(SpectralKit.grid(sp.basis)) | |
| end | |
| function residuals(coefs) | |
| # this is a mutable struct, overwrite coefficients | |
| sp.coefs = coefs | |
| c = costate(sp) | |
| 𝑑c = 𝑑costate(sp) | |
| kdot(k) = f(k, rp.α) - rp.δ * k - c(k) | |
| cdot_approx(k) = 𝑑c(k) * kdot(k) | |
| distance(sp, cdot_approx, cdot) | |
| end | |
| result = optimize(residuals, zeros(length(sp.coefs)), BFGS()) | |
| sp.coefs = result.minimizer | |
| c = costate(sp) | |
| k_range = range(0.25*kss, stop=1.75*kss, length=100) | |
| c_range = c.(k_range) | |
| plot(k_range, c_range, label="Costate", xlabel="Capital (k)", ylabel="Costate (c)", title="Saddle Path Costate") |
Generally, I think that the approach is correct: you map coefficients to residuals, and minimize that. Some technical remarks:
basis = Chebyshev(EndpointGrid(), n) ∘ BoundedLinear(k_min, k_max)
SpectralKit.jl has infinite domain bases, SemiInfRational, just center that around the steady state and you don't need to pick a min or a max. But then use InteriorGrid().
Is this discrete or continuous time? Looks like something in between.
You can implement distance as something like
mean(x -> abs2(f1(x) - f2(x)), SpectralKit.grid(sp.basis))(using Statistics, which has mean), so it is type stable. This matters if you want automatic differentiation in optimize, which you will want for larger models.
Awesome, thanks! Now I am playing around with different grids. It is the best way to learn about the method.
You may want to calculate the residuals on a denser grid when optimizing, to minimize the Runge phenomenon.
@tpapp