Created
August 21, 2025 15:43
-
-
Save korenmiklos/3c1b8cdc19fcf27b4da80631882c59e6 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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") |
Author
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.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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 useInteriorGrid().Is this discrete or continuous time? Looks like something in between.
You can implement
distanceas something like(
using Statistics, which hasmean), so it is type stable. This matters if you want automatic differentiation inoptimize, which you will want for larger models.