HW 4: Replication of Rust (1987)

Hasan K. Tosun, December 2017

In this homework, we are asked to estimate the model described in Rust (1987). However, by using a different method than the one in the original paper.

I'll use the method described in the paper first. Then, I'll move on to the "stata" estimator. (I don't know why they name it this way.)

Before going into detail, let me sketch what I'm going to do.

I'll use the data provided by Prof. Amil Petrin.

In the first section, I'll do the following:

  1. Pick a functional form for the cost function.
  2. Estimate $p(x'\mid x,d)$ by using the data.
  3. Find the fixed point of the contraction for $EV_\theta$.
  4. Find $P(d\mid x,\theta)$ by using $EV_\theta$.
  5. Calculate the likelihood function.
  6. Find parameter values $\theta$ that maximizes log-likelihood function.

In the second section, I'll follow these steps:

  1. Estimate $p(x'\mid x,d)$ by using the data.
  2. Estimate $p(d\mid x)$.
  3. Iterate on the fixed point relation to get $v(x,1)$.
  4. Use $v(x,1)$ to get $v(x,d)$.
  5. Find $u(x,d)$ by using $v(x,d)$.

Let's get started.

Part 1: Estimation of Parameters by Maximizing Partial Loglikelihood

Step 0: Import the necessary packages, import the data, set some parameters.

In [2]:
import numpy as np
import pandas as pd
from scipy.optimize import minimize
from numba import jit
import time

data = pd.read_csv("rust.csv")

cts = data['frag'].value_counts(sort = False)
nrow = data.shape[0]
pi = cts/nrow

beta = 0.9999

n = 90
x = np.array(range(1,n+1))
x.shape = (n,1)

tol = 1e-6

Step 1: Define a function that will calculate the cost of operating in state $x$, given the structural parameter $\theta_1$. I pick the linear formulation as the functional form for the cost.

In [3]:
def calc_cost(theta1):
    c = theta1 * x

Step 2: Define a function that will construct the probability transition matrix $P$. This is an $n\times n$ matrix, governing the motion of $x$'s.

In [4]:
def probs():
    P = np.zeros((n,n))
    for i in range(n-2):
        P[i,i] = pi[0]
        P[i,i+1] = pi[1]
        P[i,i+2] = pi[2]
    P[n-2,n-2] = pi[0]
    P[n-2,n-1] = 1 - pi[0]
    P[n-1,n-1] = 1

Step 3: Define the contraction function, given the structural parameters $\theta_1$ and $tr$, total cost for replacement.

In [5]:
def inside_contract(theta1,tr,ev):
    c = calc_cost(theta1)
    temp = np.exp(-c + beta * ev - ev) + np.exp(-tr - c[0,0] + beta * ev[0,0] - ev)
    return (np.log(temp) + ev)

def contraction(theta1, tr):
    dif = 1
    ev = np.zeros((n,1))
    while (dif > tol):
        ev1 = P @ inside_contract(theta1,tr,ev)
        dif = (abs(ev1-ev)).max()
        ev = ev1

Step 4: Calculate $P(0\mid x; \theta_1,tr)$.

In [6]:
def P0x(theta1, tr):
    c = calc_cost(theta1)
    ev = contraction(theta1, tr)
    pk = 1/(1 + np.exp(c - beta * ev - tr - c[0,0] + beta * ev[0,0]))

Step 5: Calculate the partial loglikelihood function: $$ \begin{equation} LL(d,x;\theta_1,tr) = \sum_{m=1}^M \sum_{t=1}^T \log P(d_{mt}\mid x_{mt};\theta_1,tr) \end{equation} $$

In [7]:
def part_likelihood(params):
    theta = params[0]
    tr = params[1]
    #print("{}  {}".format(theta,tr))
    temp = 0
    p0x = P0x(theta, tr)
    for j in range(8260):
        x_ind = data.bxt[j] - 1
        i = data.i[j]
        temp = temp + np.log(p0x[x_ind]*(i==0) + (1-p0x[x_ind])*(i==1)+0.0000001)
    temp = temp[0,]

Step 6: Maximization of Partial Loglikelihood Function

$$ \begin{equation} \max_{\theta_1,tr} LL(d,x;\theta_1,tr) \end{equation} $$
In [8]:
P = probs()
init = np.array([0.001,10])

t = time.time()
bnds = ((0.0001, 0.1), (1, 40))
soln = minimize(part_likelihood, init, method="L-BFGS-B", bounds = bnds, options={'disp': True})
elapsed_time = time.time()-t
In [9]:
print("Time elapsed:\t",elapsed_time)
print("Loglikelihood:\t", -soln.fun)
print("theta1:\t\t", soln.x[0])
Time elapsed:	 445.2776126861572
Loglikelihood:	 -299.18982815
theta1:		 0.0026573741074
TR:		 9.80632865721