GARCH(1,1) MLE

Question

import pyomo.environ as pyo
import numpy as np


def get_model_garch11mle(vec_retn: np.array) -> pyo.ConcreteModel:
    mdl = pyo.ConcreteModel(name='GARCH11')
    mdl.alpha_0 = pyo.Var(bounds=(1e-6,1), initialize=.4)
    mdl.alpha_1 = pyo.Var(bounds=(1e-6,1), initialize=.4)
    mdl.beta_1 = pyo.Var(bounds=(1e-6,1), initialize=.4)
    mdl.T = range(len(vec_retn))
    mdl.h_t = pyo.Expression(mdl.T)
    for i in mdl.T:
        if i == 0:
            mdl.h_t[i] = pyo.Expression(expr=mdl.alpha_0 / (1 - mdl.alpha_1 - mdl.beta_1))
        else:
            mdl.h_t[i] = pyo.Expression(expr=mdl.alpha_0 + mdl.alpha_1 * (vec_retn[i - 1]**2) + mdl.beta_1 * mdl.h_t[i - 1])
    
    mdl.obj = pyo.Objective(expr=sum([pyo.log(mdl.h_t[i]) + vec_retn[i]**2/mdl.h_t[i] for i in mdl.T]))
    return mdl

vec_retn = np.random.randn(100)**2 /3 + np.random.randn(100)**3 /6
temp_mdl = get_model_garch11mle(vec_retn)
opt = pyo.SolverFactory('gurobi')
temp_res = opt.solve(temp_mdl)

在这个片段中，我需要使用变量(a0、a1、b1)和ht-1来更新ht的seq。好奇的是，为Pyomo写这个问题的正确方法是什么？目前，我得到了：

ValueError: Accessing the expression of expression 'ScalarExpression' 
before the Expression has been constructed (there is currently no value to return).

当删除冗余的“表达式()”时，如：

import pyomo.environ as pyo
import numpy as np


def get_model_garch11mle(vec_retn: np.array) -> pyo.ConcreteModel:
    mdl = pyo.ConcreteModel(name='GARCH11')
    mdl.alpha_0 = pyo.Var(bounds=(1e-6, 1), initialize=.4)
    mdl.alpha_1 = pyo.Var(bounds=(1e-6, 1), initialize=.4)
    mdl.beta_1 = pyo.Var(bounds=(1e-6, 1), initialize=.4)
    mdl.T = range(len(vec_retn))
    mdl.h_t = pyo.Expression(mdl.T)
    for i in mdl.T:
        if i == 0:
            mdl.h_t[i] = mdl.alpha_0 / (1 - mdl.alpha_1 - mdl.beta_1)
        else:
            mdl.h_t[i] = mdl.alpha_0 + mdl.alpha_1 * vec_retn[
                i - 1]**2 + mdl.beta_1 * mdl.h_t[i - 1]

    mdl.obj = pyo.Objective(expr=sum(
        pyo.log(mdl.h_t[i]) + vec_retn[i]**2 / mdl.h_t[i] for i in mdl.T))
    return mdl


vec_retn = np.random.randn(100)**2 / 3 + np.random.randn(100)**3 / 6
temp_mdl = get_model_garch11mle(vec_retn)
opt = pyo.SolverFactory('mosek')
temp_res = opt.solve(temp_mdl)

使用“mosek”时：

DegreeError: MOSEK does not support expressions of degree None.

使用“gurobi”时：

RuntimeError: Cannot write legal LP file.  Objective 'obj' has nonlinear terms that are not quadratic.

Answer 1

问题在于for-循环中，您不需要声明新的Expression组件，只需将表达式分配给适当的索引：

for i in mdl.T:
    if i == 0:
        mdl.h_t[i] = mdl.alpha_0 / (1 - mdl.alpha_1 - mdl.beta_1)
    else:
        mdl.h_t[i] = mdl.alpha_0 + mdl.alpha_1 * (vec_retn[i - 1]**2) + mdl.beta_1 * mdl.h_t[i - 1]