I am trying to use JuMP to solve a dynamic programming problem, but I cannot get MPEC (with Ipopt) to work.

I am solving a prediction problem. There are 2 products, “safe” and “risky”. The safe product always yields utility 0 and the risky product yields utility `β x + ξ`

, where `x`

is a (time-varying) observable characteristic, `β`

a regression coefficient and `ξ`

a (time-varying) utility shock. In each period, if the decision maker opts for the safe product, they get utility 0 and exit the market. If the decision maker opts for the risky one, then with probability `q`

(time-varying) they get it and exit the market, and with probability `1-q`

they fail and advance to the next period.

I want to minimize the sum of squared utility shocks `ξ`

, subject to the constraint that the model-implied market shares `s(β, ξ)`

match the data market shares `s0`

. Here is a minimal (non-working) example:

```
using JuMP, Ipopt, Random
# generate fake data
N = 100 # number of periods
β = 1.0 # fundamental parameter
x = randn(N) # exog. regressor
ξ0 = 0.25 * randn(N) # good-specific iid shock; in fact E[ξ] ≜ 0
u = x * β .+ ξ0 # mean utility of good 1
q = rand(N) # probability of winning when choosing good 1 (good 2 always has q = 1, u = 0)
v, V = zeros(N), zeros(N) # initialize choice-specific value function (1) and value function
for i=1:N
v = q .* u + (1 .- q) .* vcat(V[2:end], 0.0) # solve cvf(1) for all t
V = log.(1 .+ exp.(v)) # solve emax for all t
end
s0 = exp.(v) ./(1 .+ exp.(v)) # generate observed market shares
# code up MPEC
function obj(ξ...)
sum(ξ[i] .^ 2 for i=1:N)
end
function getMeanUtilities(β, ξ::Vararg{TT, NN} where NN) where {TT <: Union{Real, VariableRef}}
# recover mean utilities
u = Vector{Any}(undef, N)
for i=1:N
u[i] = x[i] * β + ξ[i]
end
u
end
function getShares(u::Vararg{TT, NN} where NN) where {TT <: Union{Real, VariableRef}}
# # initialize and recover cvf
v = Vector{TT}(undef, N)
V = Vector{TT}(undef, N)
s = Vector{TT}(undef, N)
for i=1:N
v[i] = q[i] * u[i]
V[i] = log(1 + exp(v[i]))
end
for t=1:N
for i=1:N
v[i] = q[i] * u[i] + (1-q[i]) * (i == N ? 0 : V[i+1])
V[i] = log(1 + exp(v[i]))
end
end
for i=1:N
s[i] = exp(v[i])/(1 + exp(v[i]))
end
s
end
# now define model
mdl = Model(Ipopt.Optimizer)
@variables(mdl, begin
β
ξ[1:N]
u[1:N]
end
)
register(mdl, :obj, N, obj; autodiff=true)
register(mdl, :getMeanUtilities, N+1, getMeanUtilities;
autodiff=true)
register(mdl, :getShares, N, getShares; autodiff=true)
@constraint(mdl, [i=1:N], getMeanUtilities(β, ξ...)[i] == u[i])
@NLconstraint(mdl, [i=1:N], getShares(u...)[i] == s0[i])
```

The error I get when trying to define the nonlinear constraint is the following:

```
MethodError: Cannot `convert` an object of type Array{GenericAffExpr{Float64,VariableRef},1} to an object of type VariableRef
```

How should I proceed?