Homework on estimating dynamic discrete choice models
02 September, 2021
In this homework we want to estimate a simple dynamic discrete choice model by Keane and Wolpin that we saw in class. They propose a simple model of female labor supply as follows. Utility for married woman \(i\) in period \(t\) from working (option 1) vs not working (option 0) with \(n_i\) small children is
\[ \begin{align} U_{it}^1 &= y_{it} + w_{it} - \pi n_{it} \\ U_{it}^0 &= y_{it} + x_{it}\beta + \epsilon_{it} \end{align} \] where \(y_{it}\) is the husband’s income, and importantly, \(n\subset x\), i.e. you get utility from children while not working. Let’s write the difference in those utilities as \(U_{it}^1 - U_{it}^0\),
\[ v_{it}(x_{it},w_{it},n_{it},\epsilon_{it}) = w_{it} - \pi n_{it} -x_{it}\beta - \epsilon_{it} \] and define the work indicator \(d_{it} = \mathbf{1}[U_{it}^1 > U_{it}^0]\). Next, we define the state space as observed by the econometrician as \(\Omega_{it} = (x_{it},w_{it},n_{it})\), i.e. we don’t observe \(\epsilon\), but the decision maker does.
Woman \(i\) will work in \(t\) if \(U_{it}^1 > U_{it}^0\), i.e. if \(v_{it}(x_{it},w_{it},n_{it},\epsilon_{it}) > 0\), and at \(v_{it}(x_{it},w_{it},n_{it},\epsilon_{it}) = 0\) she is indifferent. Call the \(\epsilon\) that solves this the critical epsilon \(\epsilon^*(\Omega_{it})\). We have that, given state \(\Omega_{it}\),
\[ i\text{ chooses to }\begin{cases}\text{work in t}&\text{if }\epsilon_{it} < \epsilon^*(\Omega_{it}) \Rightarrow U_{it}^1 > U_{it}^0 \\ \text{not work in t}& \text{if }\epsilon_{it} > \epsilon^*(\Omega_{it}) \Rightarrow U_{it}^1 < U_{it}^0\end{cases} \]
Assuming that \(\epsilon\) is independent of \(\Omega\), the probability of working is computed by looking at each \(i\)’s \(\epsilon_{it}\) and couting whether it’s below \(\epsilon^*(\Omega_{it})\):
\[ \Pr[d_{it}=1|\Omega{it}] = \int_{-\infty}^{\epsilon_{it}^*} dF_{\epsilon}(\epsilon|\Omega{it}) = \int_{-\infty}^{\epsilon_{it}^*} dF_{\epsilon}(\epsilon) \]
Clearly \(\Pr[d_{it}=0|\Omega{it}] = 1- \Pr[d_{it}=1|\Omega{it}]\), and so the likelihood for a random sample of \(N\) females observed for \(T\) periods is
\[ L(\beta,\pi,F_{\epsilon};x_{it},w_{it}) = \Pi_{i=1}^N \Pi_{t=1}^T \Pr[d_{it}=1|\Omega{it}]^{d_{it}} \Pr[d_{it}=0|\Omega{it}]^{1-d_{it}} \]
Simulating a data set
We start out with a simulated data set that we will use in estimation.
# true parameters
p = list()
p$beta1 = 1
p$beta2 = 0.5
p$beta3 = 0.5
p$pi = 1.5
p$N = 2000
p$T = 60 # final year of observation
get.data <- function(p){
d = data.table(expand.grid(id = 1:p$N,it=1:p$T))
d[,n := sample(c(0,1,2),size=1), by=id]
d[,x1 := runif(n=nrow(d))]
d[,x2 := runif(n=nrow(d))]
d[,w := rnorm(nrow(d),mean=5)] # suppose we observe wages of non-workers
d[,eps := rnorm(mean=0,sd=1,n=nrow(d))]
d[,v:= w - p$pi*n - p$beta1*x1- p$beta2*x2 - p$beta3*n -eps]
d[,work := v > 0]
flog.info("participation rate: %f",d[,mean(work)])
return(d)
}
d <- get.data(p)## INFO [2021-09-02 21:49:19] participation rate: 0.826733R-golfing
Question 1. Write a one-liner that computes the means of x,w,eps,v and work from this data.table. count how many characters you had to type with nchar("YOUR_CODE_IN_QUOTES_LIKE_THIS"). I have 32.
nchar("d[,lapply(.SD,mean),.SDcols=4:9]")## [1] 32d[,lapply(.SD,mean),.SDcols=4:9]## x1 x2 w eps v work
## 1: 0.5007173 0.5001701 4.991227 0.001809978 2.175614 0.8267333Estimating a static model
Question 2. Identification.
- Run a regression of
vonw,x,nand see if you can recover the true parameters. - Make suitable assumptions to be able to recover all parameters.
summary(lm(v~w+x1+x2+n,d))##
## Call:
## lm(formula = v ~ w + x1 + x2 + n, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.2727 -0.6818 -0.0011 0.6750 4.1138
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.013447 0.016740 -0.803 0.422
## w 1.005841 0.002889 348.218 <2e-16 ***
## x1 -1.004006 0.010022 -100.179 <2e-16 ***
## x2 -0.527031 0.010019 -52.604 <2e-16 ***
## n -2.001931 0.003544 -564.896 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.002 on 119995 degrees of freedom
## Multiple R-squared: 0.7899, Adjusted R-squared: 0.7899
## F-statistic: 1.128e+05 on 4 and 119995 DF, p-value: < 2.2e-16p$pi = 0
d = get.data(p)## INFO [2021-09-02 21:49:20] participation rate: 0.993817summary(lm(v~w+x1+x2+n,d))##
## Call:
## lm(formula = v ~ w + x1 + x2 + n, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0506 -0.6696 -0.0004 0.6699 4.3195
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.002358 0.016622 0.142 0.887
## w 0.996765 0.002872 347.073 <2e-16 ***
## x1 -0.976390 0.009983 -97.803 <2e-16 ***
## x2 -0.499000 0.009978 -50.008 <2e-16 ***
## n -0.498897 0.003548 -140.623 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9961 on 119995 degrees of freedom
## Multiple R-squared: 0.5587, Adjusted R-squared: 0.5587
## F-statistic: 3.797e+04 on 4 and 119995 DF, p-value: < 2.2e-16Question 3. Functional Forms.
- What does the likelihood function above remind you of?
- Sometimes people use a linear probability model to approximate a probit model. Which parameters can we recover using an LPM here? How well does it predict the average probability of work? Give the quartiles of the predicted probabilities.
- Then do the same exercise with a probit model.
lmod = lm(work ~w+ x1 + x2 +n,d)
d[,lpm := predict(lmod,d) ]
summary(lmod)##
## Call:
## lm(formula = work ~ w + x1 + x2 + n, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99374 -0.00286 0.00599 0.01476 0.05472
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9499237 0.0012900 736.399 < 2e-16 ***
## w 0.0116735 0.0002229 52.376 < 2e-16 ***
## x1 -0.0117696 0.0007748 -15.191 < 2e-16 ***
## x2 -0.0060548 0.0007744 -7.819 5.37e-15 ***
## n -0.0053893 0.0002753 -19.574 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0773 on 119995 degrees of freedom
## Multiple R-squared: 0.02763, Adjusted R-squared: 0.02759
## F-statistic: 852.3 on 4 and 119995 DF, p-value: < 2.2e-16d[,list(mean(work),mean(lpm))]## V1 V2
## 1: 0.9938167 0.9938167d[,quantile(lpm)]## 0% 25% 50% 75% 100%
## 0.9387564 0.9849863 0.9938166 1.0026085 1.0469731probit model:
gmod = glm(work ~ w+ x1 + x2 +n,family =binomial(link="probit"),d)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredd[,glm := predict(gmod,d,type="response")]
d[,list(mean(work),mean(lpm),mean(glm))]## V1 V2 V3
## 1: 0.9938167 0.9938167 0.9938075summary(gmod)##
## Call:
## glm(formula = work ~ w + x1 + x2 + n, family = binomial(link = "probit"),
## data = d)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.3668 0.0020 0.0114 0.0479 1.5267
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.08476 0.09726 0.871 0.384
## w 1.01525 0.02505 40.534 <2e-16 ***
## x1 -1.06224 0.06805 -15.609 <2e-16 ***
## x2 -0.55302 0.06389 -8.655 <2e-16 ***
## n -0.52773 0.02753 -19.169 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 9026.9 on 119999 degrees of freedom
## Residual deviance: 5426.5 on 119995 degrees of freedom
## AIC: 5436.5
##
## Number of Fisher Scoring iterations: 10d[,quantile(glm)]## 0% 25% 50% 75% 100%
## 0.1490305 0.9987606 0.9999299 0.9999978 1.0000000Let’s make the wage a bit more realistic. Suppose we knew completed years of education, and we had a variable \(z_1\) capturing potential experience (age - educ - 6), which determines the wage, \(z_2 \sim U(0,1)\) as a second explanatory variable and \(\eta\) a random error. Additionally, let’s be more reasonable and assume that for all non-workers, we don’t observe a wage, for which we code 0.
\[ w_{it} = z_{it}'\gamma + \eta_{it} \]
We posit that \(\eta\) and \(\epsilon\) are joint normal with mean zero and var-cov matrix
\[ \sum = \left[ \begin{array}{cc} \sigma_\eta^2 & \sigma_{\eta \epsilon} \\ \sigma_{\eta \epsilon}& \sigma_\epsilon^2 \end{array} \right] = \left[\begin{array}{cc} 0.5 & 0.5 \\ 0.5 & 0.8 \end{array} \right] \]
Question 4. Change the function get.data to reflect those changes. Add a new parameters \(\gamma_1=0.02,\gamma_2=0.09\) to p. Then run a wage regression and see if you can recover the \(\gamma\)s. Report the value of \(E[\eta|\text{work}]\) and of \(E[\eta]\).
p$gamma1 = 1
p$gamma2 = 2
p$var_eps = 0.5
p$var_eta = 0.8
p$cov_eps_eta = 0.5
get.data <- function(p){
set.seed(111)
d = data.table(expand.grid(id = 1:p$N,it=20:p$T))
d[,n := sample(c(0,1,2,3),size=1), by=id]
d[,x1 := runif(n=nrow(d))]
d[,x2 := runif(n=nrow(d))]
d[,educ := sample(8:18,size=1),by=id]
# d[,z1 := it - educ - 6]
# d[,z1 := runif(n=nrow(d),1,100)]
d[,z1 := runif(n=nrow(d),-5,5)]
d[,z2 := runif(n=nrow(d),-5,5)]
mv = rmvnorm(nrow(d),sigma=matrix(c(p$var_eps,p$cov_eps_eta,p$cov_eps_eta,p$var_eta),2,2))
d[,eps := mv[,1]]
d[,eta := mv[,2]]
d[,w := p$gamma1*z1 + p$gamma2*z2 + eta]
d[,v:= w - p$pi*n - p$beta1*x1- p$beta2*x2 - p$beta3*n - eps]
d[,work := v > 0]
d[!(work), w := 0]
return(d)
}
d = get.data(p)
summary(lm(w ~ z1 + z2 ,d[(work)]))##
## Call:
## lm(formula = w ~ z1 + z2, data = d[(work)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6297 -0.5972 0.0018 0.6013 3.6374
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.144675 0.010778 13.42 <2e-16 ***
## z1 0.986264 0.001950 505.81 <2e-16 ***
## z2 1.963608 0.003302 594.59 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8945 on 34749 degrees of freedom
## Multiple R-squared: 0.9258, Adjusted R-squared: 0.9258
## F-statistic: 2.166e+05 on 2 and 34749 DF, p-value: < 2.2e-16flog.info("E[eta] = %f",d[,mean(eta)])## INFO [2021-09-02 21:49:23] E[eta] = -0.000103flog.info("E[eta|work] = %f",d[work==TRUE,mean(eta)])## INFO [2021-09-02 21:49:23] E[eta|work] = 0.035422flog.info("participation rate: %f",d[,mean(work)])## INFO [2021-09-02 21:49:23] participation rate: 0.423805Question 5. Estimate the static model with maximum likelihood.
- write a likelihood function that takes our simulated dataset
d, the contents ofpand returns the value of a the likelihood function (i.e. a number)
# the way you set up this function depends a lot on the optimizer you are going to use later on. We will be using maxLik.
static.like <- function(beta1,beta2,beta3,gamma1,gamma2,var_eps,var_eta,cov,data){
# enforce non-negativity constraint on parameters.
# some optimizers allow to have box constraints on choice variables,
# which would be preferrable. Here, in the absence of such functionality,
# this is an alternative solution.
if (var_eps<0 | var_eta <0){
return(NA)
}
# compute \xi
xi = var_eps + var_eta - 2*cov
if (xi<0) {
return(NA)
}
# compute probability of no work
d = data.table(nw_index=data[,-(gamma1*z1 + gamma2*z2 - beta1*x1 - beta2*x2 - beta3*n)/sqrt(xi)])
d[,pr0 := pnorm(nw_index)]
# joint probability of work and certain wage
# pr(xi > -v and eta = w - z \gamma) =
# pr(eta-eps > -v and eta = w - z \gamma) =
# pr( eps < w - xb )
d[,pr_w := pnorm(data[,(w - beta1*x1 - beta2*x2 - beta3*n)/sqrt(var_eps)])]
# compute value of log likelihood for each individual
d[,work := data[,work]]
d[(work),ll := log(pr_w) ]
d[!(work),ll := log(pr0)]
# print(d[,mean(pr_w,na.rm=T)])
# return sum of log of likelihood
if (!is.finite(d[,sum(ll)])){
return(NA)
}
return(d[,sum(ll)])
}
# test run
x=static.like(p$beta1,p$beta2,p$beta3,p$gamma1,p$gamma2,p$var_eps,p$var_eta,p$cov,d)- check the likelihood function: produce a plot that looks like this.
# wrapper 1: a closure to include the dataset (can't pass as parameter to optimizer!)
dwrap.like <- function(beta1,beta2,beta3,gamma1,gamma2,var_eps,var_eta,cov){
static.like(beta1,beta2,beta3,gamma1,gamma2,var_eps,var_eta,cov,d)
}
# wrapper 2: maps a numeric vector to function arguments.
mwrap.like <- function(x){
static.like(beta1 =x[1],
beta2 =x[2],
beta3 =x[3],
gamma1 =x[4],
gamma2 =x[5],
var_eps =x[6],
var_eta =x[7],
cov =x[8],
data = d)
}
# now check that the likelihood is correct.
# make a grid of values for each param +- 20 % of starting value
n = 50
l = 0.2
l1 = list(beta1=p$beta1,beta2=p$beta2,beta3=p$beta3,gamma1=p$gamma1,gamma2=p$gamma2,var_eps=p$var_eps,var_eta=p$var_eta,cov_eps_eta=p$cov)
grid = data.frame(lapply(l1,function(x) seq((1-l)*x,(1+l)*x,length=n)))
df = data.frame(lapply(l1,function(x)rep(x,n)))
d0 = data.frame()
for (i in names(grid)){
tmp = data.frame(df,param=i)
tmp[,i] = grid[,i]
d0 = rbind(d0,tmp)
}
d0$LL = mapply(dwrap.like,d0$beta1,d0$beta2,d0$beta3,d0$gamma1,d0$gamma2,d0$var_eps,d0$var_eta,d0$cov_eps_eta)
m = melt(d0,measure.vars="LL")## Warning: The melt generic in data.table has been passed a data.frame and will
## attempt to redirect to the relevant reshape2 method; please note that reshape2
## is superseded and is no longer actively developed, and this redirection
## is now deprecated. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace, i.e.
## reshape2::melt(d0). In the next version, this warning will become an error.d1 = data.frame()
for (i in names(grid)){
tmp = data.frame(pname=i,pval=m[m$param==i,i],LL=m[m$param==i,]$value)
# tmp$truth = rep(p[[i]],nrow(tmp))
d1 = rbind(d1,tmp)
}
truth= data.frame(pname=names(l1),val=unlist(l1))
ggplot(d1,aes(x=pval,y=LL)) + geom_line() + geom_vline(data=truth,aes(xintercept=val),color="red")+ facet_wrap(~pname,scales = "free")
- Set up an optimizer that searches for a maximum over that function.
require(maxLik)## Loading required package: maxLik## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'maxLik'# could constrain var_eps and var_eta to be positive
# A <- matrix(c(rep(0,5),1,0,0,
# rep(0,5),0,1,0), 2, 8,byrow=T)
# B <- c(0,0)
# mle = maxLik(mwrap.like,grad=NULL,hess=NULL,start=unlist(l1) + runif(length(l1))*0.1,fixed=c("gamma2"),control=list(printLevel=1),constraints=list(ineqA=A,ineqB=B))
# prefer setting mwrap.like=NA if negative.
# mle = maxLik(mwrap.like,grad=NULL,hess=NULL,start=unlist(l1) + runif(length(l1))*0.1,fixed=c("gamma2"),control=list(printLevel=1),method="BFGS")
# summary(mle)Question 6. Do a counterfactual experiment with this model. Illustrate how labor supply varies as we introduce a child care subsidy of \(\tau=1\) dollars if the woman works while \(n_{it}>0\).
- You can use the true parameter values stored in
p. - Ideally you would show a plot with two overlaid density estimates.
- A kernel estimate of the distribution of \(\Pr(d_{it}=1|\Omega_{it},\tau=0)\)
- A kernel estimate of the distribution of \(\Pr(d_{it}=1|\Omega_{it},\tau=1)\)
- Also report a table with mean and median of that distribution by scenario.
# use the existing data in d
p$tau = 1
d[,pr1 := pnorm(v/(p$var_eps + p$var_eta - 2*p$cov_eps_eta))]
d[,pr1_tau := pnorm((v+p$tau)/(p$var_eps + p$var_eta - 2*p$cov_eps_eta))]
m = d[,list(pr_work=pr1,scenario="tau=0")]
m = rbind(m,d[,list(pr_work=pr1_tau,scenario="tau>0")])
setnames(m,c("pr_work","scenario"))
pander(m[,list(med=median(pr_work),mean=mean(pr_work)),by=scenario])| scenario | med | mean |
|---|---|---|
| tau=0 | 2.044e-07 | 0.4239 |
| tau>0 | 0.04168 | 0.4739 |
ggplot(data=m,aes(x=pr_work,color=scenario)) + geom_density()
Question 7. Add experience (\(h\)) to our dataset and to the wage equation. The law of motion of \(h\) is
\[ h_{it} = h_{it-1} + d_{it-1} \] i.e. \(h\) increases by 1 if the individual worked last period.
h_work_update <- function(w_,gamma3,v_){
# compute w = w_ + gamma3*h
# w = p$gamma1*z1 + p$gamma2*z2 + p$gamma3*h + eta
# w_ = p$gamma1*z1 + p$gamma2*z2 + eta
# compute v = w - v_
# v_ = p$pi*n - p$beta1*x1- p$beta2*x2 - p$beta3*n - eps
# compute work = v > 0
# compute h[t] = h[t-1] + work[t-1]
n = length(w_)
w = rep(0,n)
h = rep(0,n)
work = rep(FALSE,n)
# period 1: zero experience.
w[1] = w_[1] + gamma3*h[1]
work[1] = w[1] - v_[1] > 0
for (i in 2:n){
h[i] = h[i-1] + work[i-1]
w[i] = w_[i] + gamma3*h[i]
work[i] = w[i] - v_[i] > 0
}
v = w - v_
w[!work] = 0
return(list(work=work,h=h,w=w,v=v))
}
p$gamma3 = 0.9
get.data <- function(p){
set.seed(111)
d = data.table(expand.grid(id = 1:p$N,it=20:p$T))
d[,n := sample(c(0,1,2,3),size=1), by=id]
d[,x1 := runif(n=nrow(d))]
d[,x2 := runif(n=nrow(d))]
d[,educ := sample(8:18,size=1),by=id]
d[,z1 := it - educ - 6]
# d[,z1 := runif(n=nrow(d),1,100)]
d[,z1 := runif(n=nrow(d),-5,5)]
d[,z2 := runif(n=nrow(d),-5,5)]
mv = rmvnorm(nrow(d),sigma=matrix(c(p$var_eps,p$cov_eps_eta,p$cov_eps_eta,p$var_eta),2,2))
d[,eps := mv[,1]]
d[,eta := mv[,2]]
# add wage, human capital and v recursively
setkey(d,id,it) # sorts by id and it
d[,c("work","h","w","v") := h_work_update(p$gamma1*z1 + p$gamma2*z2 + eta,
p$gamma3,
p$pi*n - p$beta1*x1- p$beta2*x2 - p$beta3*n - eps),by=id]
return(d)
}
d = get.data(p)
summary(lm(w ~ z1 + z2 ,d[(work)]))##
## Call:
## lm(formula = w ~ z1 + z2, data = d[(work)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.034 -8.447 -0.290 8.141 24.042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.62925 0.03573 465.36 <2e-16 ***
## z1 0.75447 0.01234 61.16 <2e-16 ***
## z2 1.48808 0.01257 118.37 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.691 on 74650 degrees of freedom
## Multiple R-squared: 0.1864, Adjusted R-squared: 0.1864
## F-statistic: 8552 on 2 and 74650 DF, p-value: < 2.2e-16flog.info("E[eta] = %f",d[,mean(eta)])## INFO [2021-09-02 21:49:41] E[eta] = -0.000103flog.info("E[eta|work] = %f",d[work==TRUE,mean(eta)])## INFO [2021-09-02 21:49:41] E[eta|work] = 0.023550flog.info("participation rate: %f",d[,mean(work)])## INFO [2021-09-02 21:49:41] participation rate: 0.910402Question 8. Check the created dataset with at least those two conditions:
- everybody starts with 0 experience
work = v > 0everywhere
library(testthat)
test_that("everybody starts with 0 experience", {
expect_true(all(d[it==min(it),h==0]))
})## Test passed 🌈test_that("work = v > 0", {
expect_true(all(d[,work == (v > 0)]))
})## Test passed 🥳Estimating a Dynamic Model
Now that \(h\) has been added, our model has finally become dynamic: Tomorrow’s value depends on what you do today. Here, tomorrow’s wage depends on whether you worked in all previous periods, hence trading off the decision of whether or not to work has a component that describes the future consequences of each choice.
Question 9. solve the dynamic model.
- write down the deterministic state space of the value function, \(\Omega^-\), i.e. the one without unobservalbe shocks.
- think about how to represent that. (\((x1,x2,gamma1,gamma2)\) are all continuous variables!)
- Realize that 2. is a problem. We need to evaluate the (difference in) value functions at each point in the dataset. so we would need to interpolate our value function many times.
- doing 3. requires us to first compute the value function on some grid of points (like
x1 = [0,0.5,1]etc) and then to find out the value atx1=0.21for example. - Or we take a shortcut and just say we will approximat this large continuous space by
npoints in eachx1,x2,z1,z2. you would have to bin the data according to the categories you come up inx1,x2,z1,z2(likex1=0.21 => x1_bin_1) for example. Then you’d evaluate the likelihood function on the actual data (just like up to now), but the value function inside the likelihood on this discretized version of the state space.
p$nMC = 100
p$nxb = 50
p$nzg = 50
tens1 <- tensorFunction(R[i,j,k] ~ x[j,k] + a[i])
vfun <- function(p,d){
# time periods
nt = length(p$T:20)
# deterministic state space
d[,xb := p$beta1*x1+ p$beta2*x2 + p$beta3*n]
d[,zg := p$gamma1*z1 + p$gamma2*z2]
r = d[,list(xb=range(xb),zg=range(zg))]
determ = data.table(
expand.grid(xb=seq(from=r$xb[1],to=r$xb[2],length.out = p$nxb),
zg=seq(from=r$zg[1],to=r$zg[2],length.out = p$nzg)
)
)
# as arrays
xb = array(determ$xb,c(p$nxb,p$nzg))
zg = array(determ$zg,c(p$nxb,p$nzg))
# integration setup: monte carlo shocks for eta and eps
sigma=matrix(c(p$var_eps,p$cov_eps_eta,p$cov_eps_eta,p$var_eta),2,2)
# expected value function
EV = array(NA,c(p$nxb,p$nzg,nt,nt))
# choice specific value functions
V = array(NA,c(p$nMC,p$nxb,p$nzg,nt,nt,2))
# time loop
for (it in nt:1){
flog.info("doing period = %d",it)
# notice the relationship between age and potential experience
for (ih in 0:(it-1)){
iih = ih + 1 # get a valid (1-based) index
# we will compute this for all values of xb and zg inside a block of (ih,it) combinations
if (it<nt){
EV1 = EV[ , ,iih+1,it+1]
EV0 = EV[ , ,iih,it+1]
} else {
EV1 = array(0,c(p$nxb,p$nzg))
EV0 = array(0,c(p$nxb,p$nzg))
}
# monte carlo integration
shocks = mvtnorm::rmvnorm(p$nMC, sigma=sigma)
# print(dim(tens1(xb + EV0 ,shocks[ ,1])))
# print(iih);print(it)
# print(dim(V[ , , ,iih,it,1]))
V[ , , ,iih,it,1] = tens1(xb + EV0 ,shocks[ ,1])
V[ , , ,iih,it,2] = tens1(zg + EV1 + p$gamma3*ih,shocks[ ,2])
Vmax = pmax(V[ , , ,iih,it,1],V[ , , ,iih,it,2])
# print(dim(apply(Vmax,c(2,3),mean)))
EV[ , ,iih,it] = apply(Vmax,c(2,3),mean)
}
}
return(EV)
}
v=vfun(p,d)## INFO [2021-09-02 21:49:59] doing period = 41
## INFO [2021-09-02 21:50:13] doing period = 40
## INFO [2021-09-02 21:50:15] doing period = 39
## INFO [2021-09-02 21:50:17] doing period = 38
## INFO [2021-09-02 21:50:19] doing period = 37
## INFO [2021-09-02 21:50:21] doing period = 36
## INFO [2021-09-02 21:50:23] doing period = 35
## INFO [2021-09-02 21:50:24] doing period = 34
## INFO [2021-09-02 21:50:27] doing period = 33
## INFO [2021-09-02 21:50:28] doing period = 32
## INFO [2021-09-02 21:50:32] doing period = 31
## INFO [2021-09-02 21:50:34] doing period = 30
## INFO [2021-09-02 21:50:36] doing period = 29
## INFO [2021-09-02 21:50:37] doing period = 28
## INFO [2021-09-02 21:50:39] doing period = 27
## INFO [2021-09-02 21:50:40] doing period = 26
## INFO [2021-09-02 21:50:41] doing period = 25
## INFO [2021-09-02 21:50:42] doing period = 24
## INFO [2021-09-02 21:50:43] doing period = 23
## INFO [2021-09-02 21:50:44] doing period = 22
## INFO [2021-09-02 21:50:45] doing period = 21
## INFO [2021-09-02 21:50:46] doing period = 20
## INFO [2021-09-02 21:50:48] doing period = 19
## INFO [2021-09-02 21:50:49] doing period = 18
## INFO [2021-09-02 21:50:49] doing period = 17
## INFO [2021-09-02 21:50:51] doing period = 16
## INFO [2021-09-02 21:50:51] doing period = 15
## INFO [2021-09-02 21:50:52] doing period = 14
## INFO [2021-09-02 21:50:53] doing period = 13
## INFO [2021-09-02 21:50:54] doing period = 12
## INFO [2021-09-02 21:50:54] doing period = 11
## INFO [2021-09-02 21:50:55] doing period = 10
## INFO [2021-09-02 21:50:55] doing period = 9
## INFO [2021-09-02 21:50:56] doing period = 8
## INFO [2021-09-02 21:50:56] doing period = 7
## INFO [2021-09-02 21:50:56] doing period = 6
## INFO [2021-09-02 21:50:57] doing period = 5
## INFO [2021-09-02 21:50:57] doing period = 4
## INFO [2021-09-02 21:50:57] doing period = 3
## INFO [2021-09-02 21:50:57] doing period = 2
## INFO [2021-09-02 21:50:57] doing period = 1persp(v[25, , ,40])
persp(v[ , , 1,1])
Question 10. estimate dynamic model