Solving toy model with autoregressive term in delta log

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()

You can try other NLP solvers. The default NLP solver (IPOPT) does not seem to converge to global infeasibility but CONOPT does.

You can set the solver as follows:

toyOneModel.solve(solver="CONOPT")

Conopt actually gives you good information about the infeasibility of your model. Conopt deposits this in the listing file, so use the following code:

    toyOneModel.solve(solver="conopt", options=Options(equation_listing_limit=100, listing_file="toyOneModel.lst"))
    with open("toyOneModel.lst") as f: 
        print(f.read())

the relevant part of the listing file reads

 ** An equation in the pre-triangular part of the model cannot
    be solved because there are no non-fixed variables left.

    Residual=                  0.35152359
    Tolerance (Tol_Feas_Tria)= 1.00000000E-08


    Variable LA(0) fixed at Value = 0.693047

    Variable LA(1) fixed at Value = 0.693047

    Variable EA(0) fixed at Value = 0

    Equation productivity(0) cannot be solved. No free variables left.
    Residual after last transformation = 0.351524

The bounds should be clear to you because you set them explicitly, now take a look at the first equation of productivity. You get this by print(productivity.getEquationListing(n=1)) after the solve method. This gives you:

productivity(0)..  - 1.5*LA(0) + LA(1) - EA(0) =E= 0.005 ; (LHS = -0.346523587779806, INFES = 0.351523587779806 ****

and it should be clear that this is an infeasible constraint.

In order to make this all happen, you need to import Options from GAMSPy and set the options equation_listing_limit (to use getEquationListing) and listing_file to get Conopt’s messages about the infeasibility in the solve method (as the code above shows you).

Good luck!

-Michael