Design question for model with functional forms as parameters

New to Julia
1

I asked a question at stackoverflow about how to properly design a model. I’m still unsure how to do it, partly because I have this old design that no longer works and that I’m trying to rescue, and partly because I’m struggling with the basics: struct, method, functor, abstract types, concrete types. I’ve read relevant portions of the docs, but as I’m no programmer I’m still far from grasping the basics.

The way I had things set up back in Julia 0.5 times, was something like this:

  • an immutable struct to declare types for parameter values and functional forms (that depend on the parameters).
  • a method with concrete types and default values
  • This Model was then manipulated to conduct experiments to see the effect of changing parameter values and functional forms. I had these stored in separate files and would load them in a loop and produce tables of results for each (each combination of parameters and functional forms corresponded to some reasonably well-known paper in the field).
  • I was able to write something like sol = solve(Model(p = p1)), where p1 was one particular combination of parameters and functional forms.

That code broke and I’m in the process of rewriting it. I would like to adopt the most appropriate design. But after discussing the issue over at stackoverflow, it seems that fixing the old design is difficult and possibly not advisable. Below is where I’m currently at: a model with one parameter and one function, which I want to be able to change (idiomatically and efficiently of course).

I’m not tied to this particular piece of code. What I really want is to be able to do something like Model(p = p1) or Model(p = p1, f = f1), as I don’t mind keeping my model functional forms f1 separate from my model parameters p1, if that’s recommended.

# Define an abstract type for the model's functions
abstract type ModelParameter end
abstract type ModelFunction end

# Define a struct for functional form F
Base.@kwdef struct F{ModelParameter} <: ModelFunction
    p::ModelParameter = 1.11
end
# Define a method for functional form F
(ff::F)(x) = x + ff.p  # I read a hint about using functors, but this looks ugly.

# Define a struct for functional form G
Base.@kwdef struct G{ModelParameter} <: ModelFunction
    p::ModelParameter = 2.22
end
# Define a method for functional form G
(gg::G)(x) = x - gg.p

# Wrap one function in a struct, where the function is defined inline (the design I used to have)
Base.@kwdef struct Model1{ModelParameter, ModelFunction}
   p::ModelParameter = 2.0
   f::ModelFunction = x -> x + p
end
julia> Model1(p=1.11).f(1.0)
2.11
julia> Model1(p=2.22).f(1.0)
3.22  # "WORKING!"

# Wrap one function in a struct, where the function is defined outside (as was suggested to me at stackoverflow).
Base.@kwdef struct Model2{ModelParameter, ModelFunction}
    p::ModelParameter = 2.0
    f::ModelFunction = F()
end
julia> Model2(p=1.11).f(1.0)
2.11
julia> Model2(p=2.22).f(1.0)
2.11  # NOT WORKING!

Ultimately I want to be able to manipulate the parameters and functional forms,
e.g. changing the parameter p and the function f to G instead of the default F

julia> Model2(p=2.22, f=G()).f(1.0)  # HOW CAN I GET TO DO THIS OR SIMILAR?
2

F() is initialized with p = 1.11, and you are calling the instance of F with the the input parameter 1.0, so the result is correct.

For that specific case, I think you are being redundant in the use of the functor, since you could do:

julia> Base.@kwdef struct Model3
           p = 2.0
       end
Model3

julia> (f::Model3)(x) = x + f.p

julia> m3 = Model3()
Model3(2.0)

julia> m3(1.0)
3.0

julia> m3(2.0)
4.0

julia> m4 = Model3(-1.0)
Model3(-1.0)

julia> m4(1.0)
0.0

julia> m4(2.0)
1.

You could use:

julia> Base.@kwdef struct Model2{ModelParameter, ModelFunction}
           p::ModelParameter = 2.0
           f::ModelFunction = F(p)
       end

julia>  Model2(p=1.11).f(1.0)
2.1100000000000003

julia> Model2(p=2.22).f(1.0)
3.22

(but p is redundant)

Or, more simply perhaps:

abstract type ModelFunction end

Base.@kwdef struct F{T} <: ModelFunction
  p::T = 1.11
end
(f::F)(x) = f.p - x


Base.@kwdef struct G{T} <: ModelFunction
  p::T = 1.11
end
(g::G)(x) = g.p + x

struct Model2{T<:ModelFunction}
  f::T
end


julia> Model2(F(1.0)).f(1.0)
0.0

julia> Model2(F(2.0)).f(1.0)
1.0

julia> Model2(G(2.0)).f(1.0)
3.0

julia> Model2(G(4.0)).f(1.0)
5.0
1 Like
3

Thank you Leandro! Let me play around with the code today and get back to you. This looks promising. :smile:

P.S. When I wrote “Not Working” I meant my brain.

1 Like
4

Hi Pat,
I’ve been interested in doing something like this for econ problems, and your post sounds relevant.
Let’s see if there is overlap.

Consider a generic 2-variable utility maximization problem (UMP):
%V( \theta_{u}, \theta_{B}) = \underset{x_{1}, x_{2}}{\max} \text{ } U(x_{1}, x_{2}; \theta_{U}) \text{ s.t. } B(x_{1}, x_{2}; \theta_{B}) = 0
For convenience, stack the finite-dimensional parameters: \theta \equiv [ \theta_{U}; \theta_{B} ]
Assume: U(\cdot), B(\cdot), \theta \in \Theta satisfy regularity conditions for a unique interior solution.
Solutions: x_{1}^*(\theta), x_{2}^*(\theta), V(\theta) = U(x^*; \theta_{U})
Comparative Statics: \frac{\partial x_{i}^* }{ \partial \theta_j} (Implicit Function Theorem) and \frac{\partial V }{ \partial \theta_j} (Envelope Theorem)
This boils down to solving a system of nonlinear equations (symbolically or numerically):
F(x;\theta, U, B) \equiv \begin{pmatrix} \frac{ U_{1} }{ U_{2} } - \frac{ B_{1} }{ B_{2} } \\ B \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} where:

  • x \in X \subset \mathbb{R}^{N_x} (in this example N_x =2)
  • \theta \in \Theta \subset \mathbb{R}^{N_U + N_B}
  • U \in \mathbb{R}^{X \times \Theta_{U} } (B^A denotes the set of functions from A to B)
  • B \in \mathbb{R}^{X \times \Theta_{B} }
  • F: X \times \Theta \times \mathbb{R}^{X \times \Theta_{U} } \times \mathbb{R}^{X \times \Theta_{B} } \to \mathbb{R}^{N_x}

Challenge: F depends on finite dimensional parameters \theta and on infinite dimensional parameters U(\cdot), B(\cdot) which themselves depend on \theta.
My goal: to write a julia program/pkg that takes U(\cdot), B(\cdot), \Theta as inputs, and produces x^*(\theta),V(\theta) and comparative statics as outputs.

Example:
U(x_{1}, x_{2}; \alpha)= \alpha \log(x_1) + (1-\alpha)\log(x_2)
B(x_{1}, x_{2}; p_1, p_2, I) = I - \left( p_1 x_1 + p_2 x_2 \right)
\theta =\begin{pmatrix} \alpha , p_1 , p_2 , I \end{pmatrix}' \in \Theta \equiv (0,1)\times \mathbb{R}_{++}\times \mathbb{R}_{++}\times \mathbb{R}_{++}
Sol: x_{1}^* = \frac{\alpha I }{p_{1} }, x_{2}^* = \frac{(1-\alpha) I }{p_{2} }, V = \alpha \log\left( \frac{\alpha I }{p_{1} } \right) + (1-\alpha)\log\left(\frac{(1-\alpha) I }{p_{2} }\right)

Given the very large variety of functional forms for utility functions/budget constraints used in econ, being able to solve these types of problems generically could be very valuable.

Is this the type of problem (finite & infinite dim parameters) you’re trying to solve?

1 Like
5

Hi Albert. Yes, pretty much. I haven’t done symbolics with Julia, but I have solved problems like the ones you describe using Maple, with which it is very easy to do. I’d love to be able to do that with Julia indeed. As an aside, one of the interesting things about programming problems like this, I think, is to explore non-convex problems with multiple solutions (solving the well-behaved problems by hand with pencil and paper is pretty much covered in standard textbooks, but finding weird combos of preferences and constraints is exciting, especially for students). For numerical solutions, I think your problem and my problem fall within a very general category. I’m still learning the basics, but it would seem that there ought to be very similar optimal control problems already packaged. My interest is more about dynamics and distributions (a dynamic version of the problem you describe), so I’m looking at the DifferentialEquations package (my code was written before that, so my current plan is to refactor so that I can use that package where needed instead of the discrete routines I had naively adapted from Matlab). One of the things I want to do is fit multi-parameter distribution functions to a problem with multiple agents, hence my interest in learning how to feed different functional forms to the problem.

You’ve seen quantecon, right? That’s a good place for economics, but some of the code was written for early Julia (and adapted from code Sargent had written for Matlab in the 1980s!), so typically each problem is specifically coded and a bit tedious to generalize once you change the utility or production functions. So yes, your project would be very useful!

1 Like
6

Hola Leandro,

So I have played around with the code you showed me and I think it will work for my purpose. I do have a couple of follow-up questions (I have moved the goal posts, as your answer solved the problem as initially stated!). In the code below, you’ll see that the same parameter appears in the function and in the model. The reason for this is something like a “response function”, so that the function parameter also affects the rest of the system elsewhere. Now the question I have is how to deal with changing the parameter “globally” to get an update of both the function, say G(), and the other stuff that happens inside Model(). I hope you will understand what I mean with the code below. It’s basically the code you wrote, but with the same parameter appearing in two places.

# Define a struct
abstract type ModelFunction end

Base.@kwdef struct F{R} <: ModelFunction
  p::R = 1.11
end
(f::F)(x) = f.p - x

Base.@kwdef struct G{R} <: ModelFunction
  p::R = 2.22
end
(g::G)(x) = g.p + x

Base.@kwdef struct Model{R <: Real, N <: Integer, F <: ModelFunction}
  p::R = 3.33
  n::N = 1  #unused, here for testing purposes
  f::F = F()
end

julia> Model().f(1.0)
0.0

julia> Model(f = F(3.33)).f(1.0)
2.33

julia> Model(f = G(3.33)).f(1.0)
4.33

julia> Model(f = G(3.33)).f(2.0)
5.33

# also works without keywords, but must pass all arguments in the correct order:
julia> Model(1.11, 1, G(3.33)).f(1.0)
4.33

# If I pass no argument to G(), the default argument is used
julia> Model(f = G()).f(1.0)
3.22

# G() is the same as G(2.22)
julia> Model(f = G(2.22)).f(1.0)
3.22

# The first p parameter has no effect: G() uses its default value
julia> Model(p = 0.0, f = G()).f(1.0)
3.22
julia> Model(p = 1.11, f = G()).f(1.0)
3.22

Is it possible to modify the code to update the parameter in Model()? So as to write Model(p = new_value, f = G()) and get the same as Model(p = new_value, f = G(new_value)).

An idea I’ve had is have two distinct parameters p1 and p2 and use @assert p1 == p2 inside Model(). Unless there’s a better way.

Also if I have several functions inside Model(), should I type them the same or should I have a specific type for each, like {..., F <: ModelFunction and G <: ModelFunction}.

Note to self: I should probably also have Model(p = 1.11, f = G(p = 2.22)).f(1.0) spit a warning like “Warning! The same parameter p is assigned different values in Model() and in G()”

1 Like
7

Let’s see if this applies to your use case …
The approach I use for solving a model that comes in hundreds of variations is basically this:

  • Each object that comes in different functional forms (e.g. utility functions) is a mutable struct that knows its parameters.
  • Each object can be initialized from a set of “switches” (e.g. fix a parameter at a certain value; switch some feature on or off).
  • A Case is a collection of those switches.
  • A Model is a collection of those objects.
  • A function initializes the Model from the Case by basically calling the initializers for all the model objects.
  • Parameters that are common to multiple objects get their own object. The initializers take care of passing those into each object. This can usually be done by passing the object that contains the parameters into the initializers (who know which parameters they actually need).

I deal with loading / saving parameters by assigning each object a unique id. A parameter is then a pair of (ObjectId, parameter name) (something like (:firm, :alpha)).
Since each object knows its id, it can pull the parameters it needs out of a list of such pairs that can be loaded / saved to a file.

The workflow is then:

  • Set the switches for all the model objects to be used.
  • The rest is automatic: initialize the Model, solve it, save the results.
3 Likes
8

Thanks Hendri, this indeed sounds perfect for my needs. Is your approach based on an extension package or is it built from the ground up? Do you know of a project of this type in an open repo so I can inspect the code and adapt it to my needs? :smiley:

9

I am not sure if you are asking it general enough to your actual problem (you have a finite number of functions, or are you willing to have an arbitrary number of function within each model)?

For a definite number of functions the redudant initialization of p is probably the most clear way to do it (I am not putting type parameters, which are important for performance, here):

julia> struct A
         p
       end

julia> struct B
         p
       end

julia> (f::A)(x) = x + f.p

julia> (f::B)(x) = x - f.p

julia> struct Model
         p
         f
         g
         Model(p,F,G) = new(p,F(p),G(p))
       end

julia> Model(1.0,A,B).g(1)
0.0

julia> Model(1.0,A,B).f(1)
2.0

If you need an arbitrary number of functions within the model, there are two possibilities:

  1. the functions will be used independently, in which case is hard to escape from an array of functions (which will be not good for performance probably). In that case my impression is that some metaprogramming would be needed to generate models with variable number of functions without type instability (I am not experienced with that).

  2. the functions will be reduced somehow (the actual objective function is the sum of the input functions, for example). Then you should use a constructor for the model which does this reduction, such as:

julia> struct A
         p
       end

julia> struct B
         p
       end

julia> (f::A)(x) = x + f.p

julia> (f::B)(x) = x - f.p

julia> struct Model2
         p
         f
         Model2(p,v::Vector) = new(p, x -> sum(F(p)(x) for F in v))
       end

julia> Model2(1.0,[ A, B ]).f(1.0)
2.0

(again, type annotations are important for performance here, but they complicate a little bit the code)

1 Like
10

Thank you Leandro!

I have two projects I’m porting to Julia. One project has a small number of known functions. So for this project it sounds like the way you show at the top should cover it. The second project has an endogenous distribution of unknown shape… so that will be a greater challenge. I’m not even sure how to deal with it. I’m hoping that putting the easier project together will teach me enough to attack the second one. One day at a time. :smiley:

I see that you’re wrapping the rhs with new: I haven’t seen that before, let me read up on it! I’ll see if I can fix my code with that!

Thanks!

You’re right, removing the type annotations makes it easier to read.

Patrick.

1 Like
11

That is necessary when one cannot differentiate the call to the default constructor from the new one from the parameters. See Constructors · The Julia Language

12

Thanks Leandro! I’ll read up on it and get back to you in a few hours, probably tomorrow. :smiley:

13

ps: To get that to be type stable might be tricky (independently if that will be useful to you or not, I was just curious because I have a project for which this might be useful). There is one way:

julia> struct A
         p::Float64
       end

julia> (f::A)(x) = x + f.p

julia> struct B
         p::Float64
       end

julia> (f::B)(x) = x - f.p

julia> struct Model2{T}
         p::Float64
         f::T
         function Model2(p,v::Vector)
           f = x -> sum(F(p)(x)::Float64 for F in v) # needed to help here by providing the output type of F
           new{typeof(f)}(p, f)
         end
       end

julia> @code_warntype Model2(1.0,[ A, B ]).f(1)
Variables
  #self#::var"#13#15"{Float64, Vector{DataType}}
  x::Int64
  #14::var"#14#16"{Int64, Float64}

Body::Float64
1 ─ %1  = Main.:(var"#14#16")::Core.Const(var"#14#16")
│   %2  = Core.typeof(x)::Core.Const(Int64)
│   %3  = Core.getfield(#self#, :p)::Float64
│   %4  = Core.typeof(%3)::Core.Const(Float64)
│   %5  = Core.apply_type(%1, %2, %4)::Core.Const(var"#14#16"{Int64, Float64})
│   %6  = Core.getfield(#self#, :p)::Float64
│         (#14 = %new(%5, x, %6))
│   %8  = #14::var"#14#16"{Int64, Float64}
│   %9  = Core.getfield(#self#, :v)::Vector{DataType}
│   %10 = Base.Generator(%8, %9)::Base.Generator{Vector{DataType}, var"#14#16"{Int64, Float64}}
│   %11 = Main.sum(%10)::Float64
└──       return %11
1 Like
14

It is built from the ground up each time, and I think it has to be.
Each object needs different parameters and switches to be constructed.
There are a couple of housekeeping packages that deal with ObjectId’s and with the (to me) surprisingly complicated task of keeping track of calibrated versus fixed parameters (and with providing those to numerical optimizers in a way those can understand).
The main package is still a private repo (until something publishable emerges). Parts of that project are split off into separate packages that are public (e.g. GitHub - hendri54/CollegeStratCollege: Colleges for the college stratification project).

15

(I have just realized you are Prof. Hendricks (and not one Hendri born in 1954!): When I was researching for my PhD I was greatly inspired by your “Taxation and Long-Run Growth” article, which at the time was still a working paper printed on paper! Your articles are a great source of inspiration!)

Yes, I can see that it is rather involved! A quick look at your public repo suggests you have put together a very extensive set of routines: This could become a super package for economists! I’ll try to draw as much inspiration from it as possible.

I will follow your general recommendations and see how far I get with my project. Right now I need to spend a few hours experimenting with Leandro’s and your tips until my model can be run again. I’ll keep you informed of my progress. After I start work on the more difficult model with endogenous distributions I’ll probably come back begging for more tips! Thank you so much for both of you.

:smiley:

16

Small world…
“Involved” is really the right term. What I do is probably way too complicated. But since it works, I have little incentive to rework it.

17

It’s taken me a while, I hope you haven’t lost interest. Please let me know if you see obvious problems and/or have suggestions. Thanks!

Below is an “architecture” that feels reasonably general for my purpose. In words, I keep simple parameters in one struct, each function and the parameters it uses together in separate structs, then I gather parameters and functions together to create a model, which in turn is sent to a solve method.

I’m not sure if nesting functions into a struct is good design. The purpose is that I want to create multiple combinations of parameters and functions without copy-pasting code. So this is what a simplified version looks like:

struct G
    a :: Float64
    b :: Float64
end
(dummy::G)(x) = x^(dummy.a) + b

struct H
    a :: Float64
    c :: Float64
end
(dummy::H)(x) = x^(-dummy.a) - c

struct P
    a :: Float64
    b :: Float64
    c :: Float64
end 
function P(; a = 1.0, b = 2.0, c = 3.0)
    P(a, b, c)
end 

struct F{G,H}
    g :: G
    h :: H
end 
function F(; parameters = P(), g = nothing, h = nothing)
    a = parameters.a
    b = parameters.b
    c = parameters.c
    g = x -> G(a,b)(x) ::Float64
    h = x -> H(a,c)(x) ::Float64
    F(g, h)
    # new{typeof(g),typeof(h)}(g,h)  # do I need new?
end 

struct Model{T}
    model :: T
    function Model(; parameters = P(), functions = F())
        model = (a = parameters.a,
                 b = parameters.b,
                 c = parameters.c,
                 g = functions.g,
                 h = functions.h)
        new{typeof(model)}(model)
    end 
end 

struct solve{T}
    sol :: T
    function solve(; model = nothing, options = "default options")
        m = model
        h = x -> m.a + m.f(x) + m.g(x) ::Float64
        sol = (h = h, model = m, options = options)
        new{typeof(sol)}(sol)
    end
end

Create a model:

julia> m0 = Model()
Model{NamedTuple{(:a, :b, :c, :g, :h), Tuple{Float64, Float64, Float64, var"#7#9"{Float64, Float64}, var"#8#10"{Float64, Float64}}}}((a = 1.0, b = 2.0, c = 3.0, g = var"#7#9"{Float64, Float64}(2.0, 1.0), h = var"#8#10"{Float64, Float64}(3.0, 1.0)))

Get a solution:

julia> s0 = solve(model = m0, options = ["options"])

Attempting to understand what monster I have created:

julia> fieldnames(Model)
(:model,)

Question: why the comma after :model ?

Any sign of type instability?

julia> code_lowered(s0.sol.h)
1-element Vector{Core.CodeInfo}:
 CodeInfo(
1 ─ %1  = Core.getfield(#self#, :m)
│   %2  = Base.getproperty(%1, :a)
│   %3  = Core.getfield(#self#, :m)
│   %4  = Base.getproperty(%3, :f)
│   %5  = (%4)(x)
│   %6  = Core.getfield(#self#, :m)
│   %7  = Base.getproperty(%6, :g)
│   %8  = (%7)(x)
│   %9  = Core.typeassert(%8, Main.Float64)
│   %10 = %2 + %5 + %9
└──       return %10
)

julia> code_typed(s0.sol.h)
1-element Vector{Any}:
 CodeInfo(
1 ─ %1 = Core.getfield(#self#, :m)::Model{NamedTuple{(:a, :b, :c, :g, :h), Tuple{Float64, Float64, Float64, var"#7#9"{Float64, Float64}, var"#8#10"{Float64, Float64}}}}
│        invoke Base.getproperty(%1::Model{NamedTuple{(:a, :b, :c, :g, :h), Tuple{Float64, Float64, Float64, var"#7#9"{Float64, Float64}, var"#8#10"{Float64, Float64}}}}, :a::Symbol)::Union{}
└──      unreachable
) => Union{}

julia> @code_warntype Model()
Variables
  #self#::Type{Model}

Body::Model{NamedTuple{(:a, :b, :c, :g, :h), Tuple{Float64, Float64, Float64, var"#7#9"{Float64, Float64}, var"#8#10"{Float64, Float64}}}}
1 ─ %1 = $(QuoteNode(var"#Model#11#12"()))
│   %2 = Main.P()::Core.Const(P(1.0, 2.0, 3.0))
│   %3 = Main.F()::Core.Const(F{var"#7#9"{Float64, Float64}, var"#8#10"{Float64, Float64}}(var"#7#9"{Float64, Float64}(2.0, 1.0), var"#8#10"{Float64, Float64}(3.0, 1.0)))
│   %4 = (%1)(%2, %3, #self#)::Core.Const(Model{NamedTuple{(:a, :b, :c, :g, :h), Tuple{Float64, Float64, Float64, var"#7#9"{Float64, Float64}, var"#8#10"{Float64, Float64}}}}((a = 1.0, b = 2.0, c = 3.0, g = var"#7#9"{Float64, Float64}(2.0, 1.0), h = var"#8#10"{Float64, Float64}(3.0, 1.0))))
└──      return %4
18

Because it is a tuple:

julia> t = (1,)
(1,)

julia> typeof(t)
Tuple{Int64}

I think I would opt for something like this as the final model:

struct Model{T,G,H}
           a::T
           b::T
           c::T
           g::G
           h::H
end
Model(; parameters = P(), functions = F()) = 
        Model(parameters.a,
              parameters.b,
              parameters.c,
              functions.g,
              functions.h)

Also, I don’t see why your solve function need to be a structure. That could only be a regular function, as it does not carry data (except that which is in Model)

19

Oh great, thanks for a super fast answer. I hadn’t realized that was the tuple notation. Oh yes I’ll organize the structure of the Model as you suggest. And you’re right solve does not need to be a struct as I will not store the solution in it. Brilliant!

2 Likes