5

I am trying to simulate a probit model using a latent variable Z of the following form:

\begin{aligned} y_{i} & = \begin{cases} 1 & \; \text{if } z_{i} > 0\\ 0 & \; \text{if } z_{i} \leq 0\\ \end{cases} \nonumber \\ z_{i} & = \boldsymbol{x}_{i} \boldsymbol{\beta} + \epsilon_{i} \\ \epsilon_{i} & \sim \mathcal{N}(0,1) \nonumber \\ \end{aligned}

And given that \begin{aligned} \beta_1 =0.7; \\ \beta_2 =-0.4;\\ x_{1} & \sim \mathcal{N}(0,2);\\ x_{2} & \sim \mathcal{N}(0,3);\\ \end{aligned}

The model I am using is based of on a modified version of the standard probit model (and is used in problems related to matching markets). However, Could someone explain the process for a standard probit model?

I have referred the following links but did not understand the concept completely:

a) https://rpubs.com/cakapourani/bayesian-binary-probit-model

b) How to do simulation of Probit link?

From link b for example, I understand the logic for 

y <- pnorm(beta0 + beta1*x1 + beta2*x2)

If I have to simulate the latent variable 'Z', how would this change (if at all. Are the two cases exactly the same?).

I am not looking for a code as a solution, just the logic of starting from the latent values will do.

Xi'an
  • 90,397
  • 9
  • 157
  • 575
Prometheus
  • 207
  • 1
  • 8

1 Answers1

5

Warning: In the referenced R code

n         <- 500
beta0     <- -1
beta1     <- 5.1
beta2     <- -0.3
x1        <- runif(n=n, min=0, max=1)
x2        <- (1:n)%%2
y         <- pnorm(beta0 + beta1*x1 + beta2*x2)

the last line

y         <- pnorm(beta0 + beta1*x1 + beta2*x2)

produces the probability $p$ that the associated latent variable $Z$ is positive, not a realisation of a Bernoulli random variable $\mathcal B(p)$. And not a realisation of $Z$.

For the former it should be

y <- rbinom(1,1,pnorm(beta0 + beta1*x1 + beta2*x2))

and for the latter

z <- rnorm(1,mean=beta0 + beta1*x1 + beta2*x2)

while both can be produced simultaneously as

z <- rnorm(1,mean=beta0 + beta1*x1 + beta2*x2)
y <- (z>0)
Xi'an
  • 90,397
  • 9
  • 157
  • 575
  • That makes sense now. Thank you once again. This was really helpful. – Prometheus Dec 27 '19 at 20:47
  • A quick follow up question. The model I am working to replicate is an extension of the standard probit model where the truncation point is not continuous. For simplicity's sake, lets say .1 and -.1 i.e. if z > 0.1 then y = 1 and if z < -.1 then 0. There are some additional calculations that ensures that z does not lie between -.1 and .1. Would z still be drawn from pnorm in that case or would something like a truncated normal distribution be more appropriate? Thanks again. – Prometheus Dec 27 '19 at 21:09