An elementary Ramsey growth model

ramsey.py

"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_ramsey.html
## LICENSETYPE: Demo
## MODELTYPE: NLP


An elementary Ramsey growth model

References:
Frank P. Ramsey, A Mathematical Theory of Saving, Economics Journal,
vol.38, No. 152, December 1928.

Erwin Kalvelagen, (2003) An elementary Ramsey growth model.
http://www.gams.com/~erwin/micro/growth.gms
"""

from __future__ import annotations

import os

import gamspy.math as gams_math
from gamspy import (
    Card,
    Container,
    Equation,
    Model,
    Ord,
    Parameter,
    Set,
    Sum,
    Variable,
)


def main():
    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
    )

    # SETS #
    t = Set(
        m,
        name="t",
        records=[f"t{t}" for t in range(1, 51)],
        description="time periods",
    )
    tfirst = Set(m, name="tfirst", domain=t, description="first interval (t0)")
    tlast = Set(m, name="tlast", domain=t, description="last intervat [T]")
    tnotlast = Set(
        m, name="tnotlast", domain=t, description="all intervals but last"
    )

    tfirst[t].where[Ord(t) == 1] = True
    tlast[t].where[Ord(t) == Card(t)] = True
    tnotlast[t] = ~tlast[t]

    # SCALARS #
    rho = Parameter(m, name="rho", records=0.04, description="discount factor")
    g = Parameter(m, name="g", records=0.03, description="labor growth rate")
    delta = Parameter(
        m,
        name="delta",
        records=0.02,
        description="capital depreciation factor",
    )
    K0 = Parameter(m, name="K0", records=3.00, description="initial capital")
    I0 = Parameter(
        m, name="I0", records=0.07, description="initial investment"
    )
    C0 = Parameter(
        m, name="C0", records=0.95, description="initial consumption"
    )
    L0 = Parameter(m, name="L0", records=1.00, description="initial labor")
    b = Parameter(
        m, name="b", records=0.25, description="Cobb Douglas coefficient"
    )
    a = Parameter(m, name="a", description="Cobb Douglas coefficient")

    # PARAMETERS #
    L = Parameter(
        m, name="L", domain=t, description="labor (production input)"
    )
    beta = Parameter(
        m,
        name="beta",
        domain=t,
        description="weight factor for future utilities",
    )
    tval = Parameter(
        m, name="tval", domain=t, description="numerical value of t"
    )

    tval[t] = Ord(t) - 1

    # The terminal weight beta(tlast) computation.
    beta[tnotlast[t]] = gams_math.power(1 + rho, -tval[t])
    beta[tlast[t]] = (1 / rho) * gams_math.power(1 + rho, 1 - tval[t])
    # display beta

    # Labor is determined using an exponential growth process.
    L[t] = gams_math.power(1 + g, tval[t]) * L0

    # Cobb-Douglas coefficient a computation.
    a = (C0 + I0) / (K0**b * L0 ** (1 - b))

    # VARIABLES #
    C = Variable(m, name="C", domain=t, description="consumption")
    Y = Variable(m, name="Y", domain=t, description="production")
    K = Variable(m, name="K", domain=t, description="capital")
    I = Variable(m, name="I", domain=t, description="investment")

    # EQUATIONS #
    production = Equation(
        m,
        name="production",
        type="regular",
        domain=t,
        description="Cobb-Douglas production function",
    )
    allocation = Equation(
        m,
        name="allocation",
        type="regular",
        domain=t,
        description="household choose between consumption and saving",
    )
    accumulation = Equation(
        m,
        name="accumulation",
        type="regular",
        domain=t,
        description="capital accumulation",
    )
    final = Equation(
        m,
        name="final",
        type="regular",
        domain=t,
        description="minimal investment in final period",
    )

    # Objective function; total utility
    utility = Sum(t, beta[t] * gams_math.log(C[t]))
    production[t] = Y[t] == a * (K[t] ** b) * (L[t] ** (1 - b))
    allocation[t] = Y[t] == C[t] + I[t]
    accumulation[tnotlast[t]] = K[t.lead(1)] == (1 - delta) * K[t] + I[t]
    final[tlast] = I[tlast] >= (g + delta) * K[tlast]

    # Bounds.
    K.lo[t] = 0.001
    C.lo[t] = 0.001

    # Initial conditions
    K.fx[tfirst] = K0
    I.fx[tfirst] = I0
    C.fx[tfirst] = C0

    ramsey = Model(
        m,
        name="ramsey",
        equations=m.getEquations(),
        problem="nlp",
        sense="MAX",
        objective=utility,
    )

    ramsey.solve()

    print("Objective Function Value:  ", round(ramsey.objective_value, 4))

    # Solution visualization
    # ----------------------

    rep = Parameter(m, name="rep", domain=[t, "*"])

    rep[t, "C[t]"] = C.l[t]
    rep[t, "Y[t]"] = Y.l[t]
    rep[t, "K[t]"] = K.l[t]
    rep[t, "I[t]"] = I.l[t]

    print("Solution:\n", rep.pivot().round(3))
    # End Ramsey


if __name__ == "__main__":
    main()