# 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 :
import numpy as np
import pandas as pd
from scipy.optimize import minimize
from numba import jit
import time

cts = data['frag'].value_counts(sort = False)
nrow = data.shape
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 :
@jit
def calc_cost(theta1):
c = theta1 * x
return(c)


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 :
@jit
def probs():
P = np.zeros((n,n))
for i in range(n-2):
P[i,i] = pi
P[i,i+1] = pi
P[i,i+2] = pi
P[n-2,n-2] = pi
P[n-2,n-1] = 1 - pi
P[n-1,n-1] = 1
return(P)


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

In :
@jit
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)

@jit
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
return(ev)


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

In :
@jit
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]))
return(pk)


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 :
@jit
def part_likelihood(params):
theta = params
tr = params
#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,]
#print("{}".format(temp))
return(-temp)


Step 6: Maximization of Partial Loglikelihood Function

$$\begin{equation} \max_{\theta_1,tr} LL(d,x;\theta_1,tr) \end{equation}$$
In :
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 :
print("Time elapsed:\t",elapsed_time)
print("Loglikelihood:\t", -soln.fun)
print("theta1:\t\t", soln.x)
print("TR:\t\t",soln.x)

Time elapsed:	 445.2776126861572
Loglikelihood:	 -299.18982815
theta1:		 0.0026573741074
TR:		 9.80632865721