Main process
Simulation for each Market
This code section conducts the simulation, given the
parameters.
Note that Since we want \(\sigma =
\sqrt{1-\rho^2}\) to be real, we transform \(\rho = \frac{1}{1+\exp(\rho^*)}\) make sure
\(\rho\) stays with \([0,1]\).
sim.process = function(A, B, delta, rho) {
# foreach market
sim = foreach(i = 1:market.num, .combine = 'rbind') %:%
#foreach draw
foreach(k = 1:sim.num, .combine = '+' ) %dopar% {
prediction = rep.int(0 , firm.n) # prediction for each draw
prediction_prev = rep.int(0 , firm.n)
U0 = rnorm(1)
U.firm = rnorm(firm.n)
n.try = 0
n.pred = 0
# keep adding up n.try until it exceed the predicted n.
while (n.try <= firm.n) {
n.pred = 0 # the num of firm under this n.try
prediction = rep.int(0 , firm.n)
# for each firm
for (i.firm in 0:(firm.n-1) ) {
X = as.matrix(dat2[i, market.regressors.i])
Z = as.matrix(dat2[i, (firm.regressors.i + i.firm)])
U = as.numeric(U.firm[1 + i.firm])
rho = 1/(1 + exp(rho))
profit = X %*% B +
Z %*% A +
U0 * rho +
U * sqrt(1 - rho ^ 2) -
delta * log(n.try)
is.enter = as.integer(profit > 0)
prediction[i.firm] = is.enter
n.pred = n.pred + is.enter
} # end for each firm
if (n.pred < n.try) {
break
}
n.try = n.try + 1
prediction_prev = prediction
} # end while
c(prediction_prev, n.try - 1)
}
sim = sim / sim.num
return(sim)
}Prediction Error
This section does the simulation process, given the vector of
parameters. Note that due to some restrictions on the optim
function, we will have to flatten the parameters into one single vector.
We transform the parameters back to a list of four(\(\alpha, \beta, \delta, \rho\)) with a
helper function split.param.
split.param = function(params, length.of.each.param){
param.list = list()
i = 1
for (l in length.of.each.param){
i.start = i
i.end = i + l - 1
new.append = params[i.start : i.end] %>% matrix(ncol = 1)
param.list = append(param.list, list(new.append))
i = i + l
}
return(param.list)
}
prediction.errors = function(params){
sim.result = do.call(sim.process, split.param(params, length.of.each.param))
real.result = dat2[, 2 : (2 + firm.n) ] # N, firm1, firm2, ...frimn
v = sim.result - real.result
return(v)
}Moment condition
Next we construct the moments, according to equation \(\ref{eq:mc}\)
moment.conditions = function(params){
v = as.matrix(prediction.errors(params)) # 7 * N
M = as.matrix(dat2[,-1])
M = cbind( rep(0, nrow(M)) ,M)
## E(V0X) , 3*1
E_v0X = v[,7] * M[,market.regressors.i] # 1* N * N * 3
E_v1Z = v[,1] * M[, c(9, firm.regressors.i + 0) ] # 1* N * N * 3
E_v2Z = v[,2] * M[, c(9, firm.regressors.i + 1) ] # 1* N * N * 3
E_v3Z = v[,3] * M[, c(9, firm.regressors.i + 2) ] # 1* N * N * 3
E_v4Z = v[,4] * M[, c(9, firm.regressors.i + 3) ] # 1* N * N * 3
E_v5Z = v[,5] * M[, c(9, firm.regressors.i + 4) ] # 1* N * N * 3
E_v6Z = v[,6] * M[, c(9, firm.regressors.i + 5) ] # 1* N * N * 3
E = cbind(E_v0X, E_v1Z, E_v2Z, E_v3Z, E_v4Z, E_v5Z, E_v6Z)
return(E)
}GMM Estimator
Finally, we construct a distance for GMM(Generalized Method of Moments). Typically, an GMM estimation is obtained by \[ \hat{\theta} = \text{argmin} \quad g(\theta)'Wg(\theta) \] Where \(g(\theta)\) is the vector of expected moments, which is calculated above, and \(W\) is a weight matrix. A robust procedure is to get \(\hat{\theta}^{(1)}\) by setting \(W = \mathbb{I}\), and then get \(\hat{\theta}^{(2)}\) by setting \(\hat{W} = \frac{1}{n} g(\hat{\theta}^{(1)})g(\hat{\theta}^{(1)})'\). However, due to computational limit, we here only provide the result for \(\hat{\theta}^{(1)}\), which is still consistant but not the most efficient result.
moment.distance = function(params, W = diag(21)){
g = params %>% moment.conditions %>% colMeans %>% as.matrix() # 21 * 1
dist =t(g) %*% W %*% g
return(dist)
}