Description of the issue
How do you define equations that include auto-regressive terms?
My equation is of the type
\Delta \log X_t = a_0 + a_1 \Delta \log X_{t-1} + \epsilon.
I wrote it in gamspy like this (LA is log of A and EA is the residual that is 0 when the equation stands exactly):
productivity[t] = LA[t.lead(1)] == prod_const + (1 + prod_ar1) * LA[t] - prod_ar1 * LA[t.lag(1)] + EA[t]
When I include such equation in a toy model (see model attached) to test it, it looks like the NLP solver does not converge to the best solution. If it were to find the solution, the objective function would be 0. Yet, it is 0.083
Here is what it returns
~/gams$ python ToyBasicOneCountry2.py
Objective Function Value: 0.083
Solution:
| t | Y[t] | K[t] | LA[t] | EY[t] | EK[t] | EA[t] |
|---|---|---|---|---|---|---|
| 0 | 98.000 | 49.000 | 0.693 | 0.000 | 0.0 | 0.000 |
| 1 | 100.000 | 49.980 | 0.693 | 0.001 | 0.0 | 0.000 |
| 2 | 102.463 | 50.982 | 0.698 | 0.000 | 0.0 | -0.201 |
| 3 | 86.168 | 52.020 | 0.505 | 0.000 | 0.0 | -0.172 |
| 4 | 67.267 | 52.435 | 0.241 | 0.000 | 0.0 | -0.115 |
GAMSPy version
GAMSPy version: 1.0.1
GAMS version: 47.6.0
gamspy_base version: 47.6.0
Code below
""" Toy One Country Model
Updated TFP equation """
# from __future__ import annotations
import gamspy.math as gams_math
import numpy as np
from gamspy.math import exp
from gamspy import Card, Container, Equation, Model, Ord, Parameter, Set, Sum, Variable
import sys
def main():
m = Container()
# Define sets related to time
t = Set(m,name='t', records=range(5), description='Time periods')
tfirst = Set(m, name="tfirst", domain=t, description="first interval (t0)")
tsecond = Set(m, name="tsecond", domain=t, description="second interval (t1)")
tforecast = Set(m,name='tforecast', domain=t, description='forecast periods')
tval = Parameter(m, name="tval", domain=t, description="numerical value of t")
tval[t] = Ord(t) - 1
tfirst[t].where[Ord(t) == 1] = True
tsecond[t].where[Ord(t) == 2] = True
tforecast[t].where[Ord(t) > 2] = True
# Equation coefficients
alpha = Parameter(m, name="alpha", records=0.33, description="elasticity of substitution of capital")
dep_cap = Parameter(m,name="dep_cap",records=0.05, description="depreciation rate of capital")
inv_cap = Parameter(m, name="inv_cap", records=0.035, description="investment to GDP ratio")
gL = Parameter(m,name="gL", records=0.02, description="labour growth rate")
prod_const = Parameter(m, name="prod_const", records = 0.005, description="productivity constant")
prod_ar1 = Parameter(m, name="prod_ar1", records = 0.5, description="productivity AR(1) coefficient")
# Initial values
K0 = Parameter(m, name="K0", records=49.00, description="capital at t=0")
K1 = Parameter(m, name="K1", records=49.98, description="capital at t=1")
LA0 = Parameter(m, name="A0", records=np.log(1.9998), description="TFP at t=0")
LA1 = Parameter(m, name="A1", records=np.log(1.9998), description="TFP at t=1")
L0 = Parameter(m, name="L0", records=49.00, description="labor at t=0")
L1 = Parameter(m, name="L1", records=49.98, description="labor at t=1")
Y0 = Parameter(m, name="Y0", records=98.00, description="GDP at t=0")
Y1 = Parameter(m, name="Y1", records=100.00, description="GDP at t=1")
# Exgoneous variables
L = Parameter(m, name="L", domain=t, description="labor (production input)")
# Labor is determined using an exponential growth process.
L[t] = gams_math.power(1 + gL, tval[t]) * L0
# Endogenous variables
K = Variable(m, name="K", domain=t, description="capital")
EK = Variable(m, name="EK", domain=t, description="residual for capital function")
Y = Variable(m, name="Y", domain=t, description="gross domestic product")
EY = Variable(m, name="EY", domain=t, description="residual for production function")
LA = Variable(m, name="LA", domain=t, description="log of total factor productivity")
EA = Variable(m, name="EA", domain=t, description="residual for productivity equation")
# EQUATIONS #
production = Equation(
m,
name="production",
type="regular",
domain=t,
description="Cobb-Douglas production function",
)
accumulation = Equation(
m,
name="accumulation",
type="regular",
domain=t,
description="capital accumulation",
)
productivity = Equation(
m,
name="productivity",
type="regular",
domain=t,
description="total factor productivity function",
)
# Objective function; total utility
sse = Sum(t, (EY[t])**2 + (EK[t])**2 + (EA[t])**2)
production[t] = Y[t] == exp(LA[t]) * (K[t] ** alpha) * (L[t] ** (1 - alpha)) * exp(EY[t])
accumulation[t] = K[t] == (1 - dep_cap + EK[t]) * K[t.lag(1)] + inv_cap * Y[t.lag(1)]
productivity[t] = LA[t.lead(1)] == prod_const + (1 + prod_ar1) * LA[t] - prod_ar1 * LA[t.lag(1)] + EA[t]
# Bounds.
K.lo[t] = 0.001
Y.lo[t] = 0.001
# A.lo[t] = 0.001
# Initial conditions
K.fx[tfirst] = K0
Y.fx[tfirst] = Y0
LA.fx[tfirst] = LA0
K.fx[tsecond] = K1
Y.fx[tsecond] = Y1
LA.fx[tsecond] = LA1
EA.fx[tfirst] = 0
EA.fx[tsecond] = 0
EA.l[t] = 0
toyOneModel = Model(
m,
name="toyOneModel",
equations=m.getEquations(),
problem="nlp",
sense="MIN",
objective=sse,
)
toyOneModel.solve()
# toyOneModel.solve(output=sys.stdout)
print("Objective Function Value: ", round(toyOneModel.objective_value, 4))
# Solution visualization
# ----------------------
rep = Parameter(m, name="rep", domain=[t, "*"])
rep[t, "Y[t]"] = Y.l[t]
rep[t, "K[t]"] = K.l[t]
rep[t, "LA[t]"] = LA.l[t]
rep[t, "EY[t]"] = EY.l[t]
rep[t, "EK[t]"] = EK.l[t]
rep[t, "EA[t]"] = EA.l[t]
print("Solution:\n", rep.pivot().round(3))
if __name__ == "__main__":
main()