Status of integration, optimization, and nl solvers

The integrals are fairly simple: \int_0^T K e^{-\rho s} ds and \int_T^\infty H(s, x) e^{-\rho s} ds. I can pick the functional form of H, but probably it will be of the form Z(s)\phi(x). I suppose these are simple enough for FastGaussianQuadrature.jl?

It is about structural estimation, but not as general as your problem. It’s not a GE model, but rather a duration model with behavioral interpretation.

That seems to call a library outside of Julia. I suppose that automatic differentiation does not work here. Also, last time I worked with Optim, it was not easily parallelizable, what about NLopt?

I’ll take a look at that, thank you.

I see that FastGaussianQuadrature integrates over a closed interval, I guess Quadrature it is then.

Yes, with (shifted) Laguerre polynomials.

You can use AD to supply the derivative, it works fine.

A lot of local algorithms are sequential and not inherently parallelizable, but you can make your objective parallelizable just fine.

1 Like

I’d like to play with that. I think it could be fastest way to integrate my functions. Can you provide a reference? Wikipedia doesn’t go in depth for the use of Laguerre polynomials in integrals.

EDIT: Found it in Judd, but if you have another reference I’d appreciate it.

If you have the Judd book I think it is the best reference, really. Miranda and Fackler has a short subsection but few details.

1 Like

GitHub - QuantEcon/Expectations.jl: Expectation operators for Distributions.jl objects is largely intended to wrap the (non-adaptive) Gaussian quadrature routines in FastGaussQuadrature.

Depending on your use case, it may be the most convenient (and all of the change-of-variables required to transform the spaces is implemented).

2 Likes

Do you use NLsolve directly or through GalacticOptim?

maybe not ready for prime time, but i have been using GitHub - JuliaNLSolvers/NLSolvers.jl: No bells and whistles foundation of Optim.jl as a replacement for both Optim (unconstrained) and NLSolve

2 Likes

It looks like your solving an econ continuous time problem.
I’ve thought about doing this in julia quite a bit.

Can you please write out your problem?

1 Like

It’s a Nash bargaining problem of time of retirement for couples (Honoré & de Paula 2018). The optimization problem is

\max _{t_{1}, t_{2}}\left(\int_{0}^{t_{1}} K_{1} e^{-\rho s} d s+\int_{t_{1}}^{\infty} H_{1}\left(s, {x}_{1}\right) D\left(s, t_{2}\right) e^{-\rho s} d s-A_{1}\right) \\ \quad \times\left(\int_{0}^{t_{2}} K_{2} e^{-\rho s} d s+\int_{t_{2}}^{\infty} H_{2}\left(s, {x}_{2}\right) D\left(s, t_{1}\right) e^{-\rho s} d s-A_{2}\right)

where D\left(s, t_{j}\right)=(\delta-1) \mathbb{1}\left(s \geq t_{j}\right)+1. The objective function simplifies to

\begin{aligned} N\left(t_{1}, t_{2}\right)=& \overbrace{\left(K_{1} \rho^{-1}\left(1-e^{-\rho t_{1}}\right)+\widetilde{H}_{1}\left(t_{1}, {x}_{i}\right)+(\delta-1) \widetilde{H}_{1}\left(\max \left\{t_{1}, t_{2}\right\}, {x}_{1}\right)-A_{1}\right)}^{\equiv I} \\ & \times \underbrace{\left(K_{2} \rho^{-1}\left(1-e^{-\rho t_{2}}\right)+\widetilde{H}_{2}\left(t_{2}, {x}_{2}\right)+(\delta-1) \widetilde{H}_{2}\left(\max \left\{t_{1}, t_{2}\right\}, {x}_{2}\right)-A_{2}\right)}_{\equiv I I} \end{aligned}

with \widetilde{H}_{i}\left(t, {x}_{i}\right)=\int_{t}^{\infty} H_{i}\left(s, {x}_{i}\right) e^{-\rho s} d s.
There are 3 cases of FOC to maximize that objective function: when t^*_1 < t^*_2, when t^*_1 > t^*_2 and when t^*_1 = t^*_2. There are known bounds on t^*_1, t^*_2 for each of those cases.

In another post I ask about how to solve those FOC with NLsolve by supplying those theoretical bounds on the optimal values t^*_1, t^*_2. I’d appreciate any suggestion you may have.

2 Likes

Actually, the solution for t_1 is independent of the solution for t_2, so I can use the (faster?) univariate methods in Roots.jl. Then I can find t^*_2 with the same univariate approach.

\tilde{H} should be an ideal candidate for quadrature with fixed nodes (just shift by t).

Personally would just write a callable that evaluates that,

struct Hquad{T,V<:AbstractVector{T}}
    nodes::V
    weights::V
    ρ::T
end

function (hq::Hquad)(t, x)
    ...
end

and unit test it (because I have a tendency to mess up simple transformations like adding t and dividing by \rho), and then it should be straightforward and work seamlessly with AD.

1 Like

Hmm, this is a very good suggestion. I had defined a function

function discountedintegral(f, r, T; nodes=5_000)
    x, w = gausslaguerre(nodes)
    exp(-r*T)/r * sum( w[i]*f(x[i]/r + T) for i in 1:nodes )
end

but that computes the nodes every time it is called. With your suggestion, it would be like this:

struct Hquad{T, V <: AbstractVector{T}}
    nodes::V
    weights::V
    r::T
end

function (hq::Hquad)(f, T)
    x = hq.nodes
    w = hq.weights
    r = hq.r
    exp(-r*T)/r * sum( w[i]*f(x[i]/r + T) for i in 1:length(x) )
end

hq = Hquad(gausslaguerre(1_000)..., 0.04)

I’d wish there was a way to prevent hardcoding the number of nodes. Perhaps with memoization?

I would just save those nodes and weights in a constant once.

And unless your H is really nasty, I think you should be fine with much less than 1000. I would start with 20; when it works, Gaussian quadrature is amazingly accurate for most purposes.

I thought the point of defining a struct was to carry those nodes and weights within the struct itself. Perhaps you mean to store the number of nodes in a const? If not, I could have saved the nodes and weights and just employ a standard function that uses those consts.

I would save a version transformed by \rho in the struct already (especially if it does not change often), but the way you propose is fine too.

If you go with Expectations.jl, you can see that this is pretty much exactly what it does. Expectations.jl/iterable.jl at master · QuantEcon/Expectations.jl · GitHub is what generates the structure. The structure is Expectations.jl/types.jl at master · QuantEcon/Expectations.jl · GitHub and pretty much the only fancy behavior is that it is overloaded to support the properties you would think of as a linear opertor (i.e. Expectations.jl/iterable.jl at master · QuantEcon/Expectations.jl · GitHub and Expectations.jl/iterable.jl at master · QuantEcon/Expectations.jl · GitHub). There is a convenience operator for “functions” in Expectations.jl/iterable.jl at master · QuantEcon/Expectations.jl · GitHub

You can see all of the operations in GitHub - QuantEcon/Expectations.jl: Expectation operators for Distributions.jl objects that it supports. I bring this up because if you are writing any custom code, it would be awesome if you instead added mini features or fixed a performance/etc. issue there instead as I suspect it does close to what you need.

2 Likes

Initially I was confused, since I was not taking an expectation (not integrating over a distribution). Now that you pointed out that there is method for regular functions, I will explore that package more in depth. Thank you.

I think the way that the current interface is setup it accepts a distribution, which of course you could pass in a uniform distribution to have it just return back the simplest uniform nodes that would be the definite integral. A more convenient interface could be written if you think it is worth it.

But a few thoughts without knowing the exact details of your problem.

  • If you are doing a discounted integral to infinity, in principle you can write it as an expectation over gauss-laguerre quadrature and an exponential distribution (look at the PDF of the exponential distribution and you will see it immediately)
    • But be warned. I have found that using the change-of-variables formula for this can lead to quadrature nodes that are too far apart close to the origin and become way too large and useless on the other side. I would take a look at the nodes that the method spits out as it would depend on your discount rate.
  • If you do that approach (or even truncate it with a large T and use Gauss-Legendre quadrature), you may want to figure out a change of variable for your \hat{H} thing so it in on the [0,\infty) domain in all cases. At that point you could reuse the same nodes and quadrature points even if you are changing the t internally to that transformed function.
  • If you doing some sort of algorithm with a fixed point and need to solving these equations already on a time grid (which you can choose out of convenience or it needs to be uniform for whatever reason) , then it is sometimes more efficient to use trapezoidal or simpsons rule. That is, you could line up the grid of what you are solving with the quadrature. But this doesn’t apply if you would need to evaluate the functions directly.
  • It is entirely possible that in this case the Expectations.jl code isn’t helpful and you can just keep around the weights/etc. in your structure as @Tamas_Papp has suggested.

As always, only worry about this stuff if you find it to be too slow in practice. If this isn’t going to be the crux of your speed issues, then it isn’t worth worrying about. But it might be here given the expressions you listed above.

1 Like