# MPEC for Starters

Think of a consumer who decides how much of a certain good to consume, given prices. The consumer's behaviour is dictated by a preference parameter $\beta$, which is unobserved. We only observe data on choices, i.e. how much was consumed, at which price of the good. We will specify an economic model, and estimate the value of the preference parameter.

The utility function is given generically as 

$$
u(c;{\color{red}\beta})
$$

and there is a structural relationship dictated by theory between
demand for consumption and price of the consumption good: demand will be such that
the marginal utility of consumption is equal to price, in other words

\begin{equation}
u_{c}(c;{\color{red}\beta})=p \quad \text{ (FOC)}
\end{equation}

Our aim here is to show how one can estimate the parameter vector
$\beta$ quite easily with MPEC. This means that we set up an estimation
procedure that is a constrainted maximization problem. Our estimate
of $\beta$ will be a parameter vector that satisfies FOC
and at the same time minimizes some criterion function.

In order to faciliate estimation, we assume that consumption is measured
with error. In that case we don't observe actual consumption, but
another variable `demand`, i.e. $q=c+\epsilon$, where $\epsilon$ follows some distribution function.
Suppose we have $N$ measurements of price-demand pairs $\left\{ q_{i},p_{i}\right\} _{i=1}^{N}$.
For the sake of this example, we assume a very simple form of the
utility function; bear in mind, however, that the actual strength
of MPEC is that one can dispense with the need to find convenient
closed form solutions, just so that the model can be solved. Our utility
will be quadratic, as in

$$
u(c;{\color{red}\beta})=c-{\color{red}\beta} c^{2}
$$

and marginal utility is 
$$ u_c = 1 - 2{\color{red}\beta} c $$

Given our assumption about measurement error, we substitute out $c$


$$ u_c = 1 - 2{\color{red}\beta} (q - \epsilon) $$

The economic model prescribes that

$$ p = u_c $$

or, in other words

$$ p_i = 1 - 2{\color{red}\beta}(q_i - \epsilon_i) $$

Our estimation will be based on minimizing the sum of squared errors, $\epsilon$
subject to complying with constraint FOC. We solve the
following problem:

$$
\min_{\left\{ \epsilon\right\} _{i=1}^{N},{\color{red}\beta}}  \sum_{i=1}^{N}\epsilon_{i}^{2}  \\
s.t.  u_{c}(q_{i}-\epsilon_{i};{\color{red}\beta})  =p_{i}
$$


We notice right away that the choice of $\beta$ does not directly influence the value of the objective function, therefore both the objective and it's gradient will be independent of $\beta$. However, every choice of $\beta$ needs to go together with a corresponding choice of $\epsilon$'s which need to be chosen s.t. the constraint is satisfied. This is not a general feature of MPEC, as in many cases the parameters of interest would appear in the objective. Also, in a more elaborate model, we would probably have some data that we would want to model to come close to; In that case our objective function would be augmented by the likelihood function of that data together with our model, or we could have a moment criterion etc.



### Task 1: write down the Lagrangian of this problem



### Task 2: What is the length of the choice vector?

* What is the length of the choice vector $x$ here? what are it's elements?
* What is the number of constraints?



### Task 3: write down the constraint function and it's gradient wrt to $x$

* Denote the constraint function $g$ as the collection of all constraints (there are $m=N$ constraints and $n=N+1$ choice variables). Write down a typical element $g_i$, i.e. what does the $i$-th constraint look like? 


* What does the gradient of $g$ look like? It's a matrix. Different solvers make different assumptions here. NLopt expects an $n$ by $m$ matrix, where `g[j,i]` should be $\frac{\partial g_i}{\partial x_j}$


### Task 4: Implement using plain vanilla NLopt

Take

* N=100 grid points
* $\beta = 0.1$ as true value
* setup this as an MPEC problem in NLopt!


```julia
using Distributions
srand(12345)
normal = Normal(0,0.01)
N = 100
price = collect(linspace(0.05,0.95,N))
beta0 = 0.1
demand0 = (1.0-price) / (2*beta0)
demand = demand0 - rand(normal,N);
```


```julia
module nlopt

    using Distributions
    using NLopt
    srand(12345)
    normal = Normal(0,0.01)
    N = 100
    price = collect(linspace(0.05,0.95,N))
    beta0 = 0.1
    demand0 = (1.0-price) / (2*beta0)
    demand = demand0 - rand(normal,N)

    # objective function
    function obj(x,g)
        # x = [e_1,...,e_N,beta] i.e. (N+1,1)
        
        # gradient = [grad w.r.t e_1,
        #            ...,grad w.r.t e_N,
        #                grad w.r.t beta] i.e. (N+1,1)
        if length(g)>0
           g[:] = vcat(2.0*x[1:(end-1)],0.0)

        end
        # value of objective
        r = sum(x[1:(end-1)].^2)
        return r
    end

    # 
    function constr(r::Vector,x::Vector,g::Matrix,n,q,p)
        if length(g) > 0
            # g has to be n by m for nlopt
            # \partial g_i / \partial x_j
            g[:,:] = cat(1,diagm(2*x[end]*ones(n)),-2*(q'.-x[1:(end-1)]'))
        end
        # value of contraints
        r[:] = 1.0 - 2*x[end].*(q.-x[1:(end-1)]) .- p  # g_i
    end
    constr_clos(r::Vector,x::Vector,g::Matrix) = constr(r::Vector,x::Vector,g::Matrix,N,demand,price)

    function run()
        opt = Opt(:LD_SLSQP,N+1)
        lower_bounds!(opt,[[-Inf for i in 1:N]...,0.0])
        upper_bounds!(opt,[[Inf for i in 1:N]...,1.0])
        min_objective!(opt,obj)
        equality_constraint!(opt,constr_clos,[1e-10 for i in 1:N])
        xtol_rel!(opt,1e-4)
        ftol_rel!(opt,1e-6)

        res = optimize(opt, vcat(rand(normal,N),0.9))
        println("beta0 = $beta0")
        println("beta  = $(res[2][end])")
        r = zeros(N)
        g =zeros(N+1,N)
        constr(r,res[2],g,N,demand,price);
        println("maximal error of constraint at solution = $(maxabs(r))")
        return res
    end

end # module
```




    nlopt



### Task 5: Implement using JuMP

* Here is where the fun starts.
* Go back to HW-constrained for some inspiration!


```julia
module jump

    using JuMP 
    using Distributions
    srand(12345)
    normal = Normal(0,0.01)
    
    function run()

        # create a model

        # define constants (N, price, etc)
        N = 100
        price = collect(linspace(0.05,0.95,N))
        beta0 = 0.1
        demand0 = (1.0-price) / (2*beta0)
        demand = demand0 - rand(normal,N)

        #define JuMP variables


        # define constraints and objective
        


        # solve
        status = solve(m)
        Dict(:obj=>getobjectivevalue(m),:beta=>getvalue(beta),:eps=>getvalue(eps))
    end
end
```

    WARNING: replacing module jump





    jump



7f3439050606e249baf8ccf7d094fefba9aa374f
ipynb
Jupyter Notebook
Notebooks/mpec-starters.ipynb
snowdj/CoursePack
f3ad25012167f580b9ea97aa7e84890986648837
MIT-0
2017-06-21T16:55:54.000Z
2021-01-04T13:44:26.000Z
Notebooks/mpec-starters.ipynb
snowdj/CoursePack
f3ad25012167f580b9ea97aa7e84890986648837
MIT-0
2020-01-23T14:13:36.000Z
2020-10-01T11:10:31.000Z
Notebooks/mpec-starters.ipynb
snowdj/CoursePack
f3ad25012167f580b9ea97aa7e84890986648837
MIT-0
2017-02-24T12:23:15.000Z
2021-11-03T07:25:38.000Z
Qwen/Qwen-72B
1. YES
2. YES