Estimate Spatial Probit

The goal of this tutorial is to show a simple reproducible example of estimating a spatial or spatio-temporal probit model in R. The code estimates the model originally derived by Beron, Murdoch, and Vijverberg (2003). This R code is based on the replication materials by Franzese, Hays, and Cook (2016), who provide similar Matlab code. The data in the demo is a straight replication of a single run of Monte Carlo analysis by Franzese, Hays, and Cook (2016). This is done for the purposes of face validity.

I prefer to use optimr, but the choice of an optimization algorithm is the prerogative of the user:

library(optimr)

Clear memory:

rm(list=ls())

First run the likelihood function. I am not going to comment various parts of this function. Please read Beron, Murdoch, and Vivenberg (2003) for technical detail.

sprobit_loglik<-function(par,y,XNT,WNT,TL,n,rand_mat) {
  r<-ncol(rand_mat)
  LL_evd<-NULL
  nt<-nrow(XNT)
  i_nt<-diag(nrow(XNT)) 
  rho<-par[(length(par)-1)]
  alpha<-par[(length(par))]
  betas<-par[1:(length(par)-2)]  
  G = (i_nt-rho*WNT-alpha*TL)  
  Z<-diag(as.vector(1-2*y))   
  vcov<-Z%*%solve(G,t(solve(G)))%*%t(Z)
  ACH<-chol(solve(vcov),pivot = TRUE)
  BCH<-solve(ACH)  
  X<-as.matrix(cbind(1,XNT))
  GX<-solve(G,(X%*%betas))
  V<--Z%*%GX  
  ln_prob<-matrix(NA, nrow=nt,ncol=r)
  for (j in seq(1,r)) {
    nu<-rep(0,nt)
    for (z in seq(1,nt)) { 
      zz<-nt-(z-1)
      if (zz==nt) {
        sumterm<-0  
      } else {
        sumterm<-BCH[zz,]%*%nu
      }
      nu0<-(1/BCH[zz,zz])%*%(V[zz,1]-sumterm)
      ln_prob[zz,j]<-log(pnorm(nu0))
      nu[zz] <-qnorm(rand_mat[zz,j])*pnorm(nu0) 
    }
  }
  LL_evd<--mean(colSums(ln_prob[(n+1):nt,]))
  cat("LL_evd:",LL_evd,fill=TRUE)
  return(LL_evd)
}

Now input your data (or run the default example). Start with inputting a row standardized n by n W matrix (in the example, it is assumed to be constant over time, but does not have to be).

W<-matrix(c(0,0,0,0,0,0,0,0.25,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0.166666666666667,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0.166666666666667,0.166666666666667,0,0,0,0,0,0,0,0,0,
0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0.166666666666667,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0.2,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0.166666666666667,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0.166666666666667,0,0,0.166666666666667,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0.2,0,0.2,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0.25,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0.166666666666667,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0.25,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.142857142857143,0.142857142857143,0,0,0,0,0,0,0,0,0,0,0.142857142857143,0,0,0,0,0,0,0,0,0,0.142857142857143,0,0,0,0,0,0,0.142857142857143,0,0,0,0.142857142857143,0,0.142857142857143,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0.25,0,0.25,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0,0,0.2,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0.25,0,0,0,0.2,0,0.2,0,0,0,0,0.2,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0.125,0,0,0,0,0,0,0,0.125,0,0.125,0.125,0.125,0,0,0,0,0,0,0,0,0,0.125,0,0,0,0,0,0,0,0,0.125,0,0,0,0,0,0.125,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0.166666666666667,0.166666666666667,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0.2,0,0.2,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0.25,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0.166666666666667,0,0.166666666666667,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0.166666666666667,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0.166666666666667,0,0,0.166666666666667,0.166666666666667,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0.125,0,0.125,0,0,0,0,0,0.125,0,0,0,0,0,0.125,0,0,0,0,0,0,0.125,0.125,0,0,0,0,0,0,0,0.125,0,0,0,0,0,0,0,0,0,0,0,0,0.125,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.333333333333333,0,0,0,0,0,0,0,0.333333333333333,0,0,0.333333333333333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0,0,0.2,0,0,0,0,0,0,0,0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0.25,0,0,0,0,0,0,0.25,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666666666667,0,0,0.166666666666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0),nrow=50,ncol=50)

Input the number of observations, here it is determined by the dimensions of W.

n<-max(nrow(W),ncol(W))

Input the number of temporal periods:

t<-21

This code converts the W matrix to its spatio-temporal version. This works, as long as your W does not vary over time. Otherwise, you may need to edit this:

i_t<-diag(t)
WNT<-kronecker(i_t,W)
i_nt1<-diag(n*(t-1))
top<- matrix(0, n,n*(t-1))
right<-matrix(0, n*t,n)
TL<-cbind((rbind(top,i_nt1)),right)

Specify an NT by k matrix of k exogenous (independent) variables. In this example, they are drawn from a uniform distribution. If you are planning on running Monte Carlo simulations, this must be done outside of the simulations (as this is the fixed part of the data).

XNT<-as.matrix(runif(n*t)*5-2)

Specify your Y. Here it is generated using a spatial probit DGP, but in normal analysis, you would just replace the following 7 lines with a line that inputs your Y, as a column (e.g y<-mydata$y).

u<-rnorm(n*t)
i_nt<-diag(n*t)
rho=0.5
phi=0.25
G<-(i_nt-rho*WNT-phi*TL)
ystar<-solve(G,(-1.5 + 3*XNT + u))
y<-ifelse(ystar>0,1,0)

This creates the importance matrix, I would leave this code as is. r is the number of draws from the importance matrix.

r<-100; 
rand_mat1<-matrix(runif(n*t*(r/2)),n*t,r/2)
rand_mat2<-1-rand_mat1
rand_mat<-cbind(rand_mat1,rand_mat2)

Input starting values:

p0<-c(-1.5, 3, 0.01, 0.01)

Now call the likelihood function within a optimization algorithm (I used -optimr-). I am also interested in timing how long it takes to get the results, so I will start the timer before the optimization and end it right after.

ptm<-proc.time()
res<-optimr(par=p0,sprobit_loglik,y=y,rand_mat=rand_mat,n=n,WNT=WNT,XNT=XNT,TL=TL,
            hessian=TRUE,control=list(maxit=100000))

## LL_evd: 398.9933
## LL_evd: 417.749
## LL_evd: 449.8205
## LL_evd: 313.1288
## LL_evd: 291.4598
## LL_evd: 227.9455
## LL_evd: 177.1651
## LL_evd: 200.0274
## LL_evd: 201.8309
## LL_evd: 326.8383
## LL_evd: 170.5147
## LL_evd: 441.2834
## LL_evd: 210.2212
## LL_evd: 246.7947
## LL_evd: 177.2863
## LL_evd: 282.1029
## LL_evd: 180.8543
## LL_evd: 164.2564
## LL_evd: 187.6951
## LL_evd: 155.4019
## LL_evd: 161.7302
## LL_evd: 181.8096
## LL_evd: 163.3172
## LL_evd: 158.3203
## LL_evd: 153.0101
## LL_evd: 161.4967
## LL_evd: 154.9096
## LL_evd: 189.1505
## LL_evd: 156.0359
## LL_evd: 149.0827
## LL_evd: 158.011
## LL_evd: 153.5835
## LL_evd: 151.1689
## LL_evd: 149.3176
## LL_evd: 149.5195
## LL_evd: 149.7188
## LL_evd: 148.0768
## LL_evd: 152.0002
## LL_evd: 148.2803
## LL_evd: 150.7515
## LL_evd: 148.7647
## LL_evd: 150.0735
## LL_evd: 148.1776
## LL_evd: 146.5058
## LL_evd: 146.7227
## LL_evd: 147.9163
## LL_evd: 146.6579
## LL_evd: 147.6896
## LL_evd: 147.289
## LL_evd: 149.3861
## LL_evd: 147.2534
## LL_evd: 147.2785
## LL_evd: 146.6868
## LL_evd: 147.3028
## LL_evd: 146.4982
## LL_evd: 147.2224
## LL_evd: 146.5381
## LL_evd: 146.3759
## LL_evd: 147.2536
## LL_evd: 146.5071
## LL_evd: 146.3513
## LL_evd: 146.0162
## LL_evd: 146.2404
## LL_evd: 147.1574
## LL_evd: 145.8428
## LL_evd: 146.2055
## LL_evd: 146.5394
## LL_evd: 146.3636
## LL_evd: 146.8267
## LL_evd: 146.3759
## LL_evd: 146.8267
## LL_evd: 145.7411
## LL_evd: 145.926
## LL_evd: 146.0329
## LL_evd: 146.7527
## LL_evd: 146.7042
## LL_evd: 145.6126
## LL_evd: 146.6828
## LL_evd: 145.8591
## LL_evd: 146.8162
## LL_evd: 146.5177
## LL_evd: 145.9793
## LL_evd: 146.3283
## LL_evd: 145.7585
## LL_evd: 146.8966
## LL_evd: 146.394
## LL_evd: 146.3931
## LL_evd: 146.2025
## LL_evd: 147.079
## LL_evd: 146.3782
## LL_evd: 146.115
## LL_evd: 146.0073
## LL_evd: 146.864
## LL_evd: 145.5339
## LL_evd: 146.1329
## LL_evd: 146.4114
## LL_evd: 146.464
## LL_evd: 146.181
## LL_evd: 146.2377
## LL_evd: 146.9365
## LL_evd: 146.3675
## LL_evd: 146.8342
## LL_evd: 146.5688
## LL_evd: 146.3305
## LL_evd: 146.9933
## LL_evd: 146.2553
## LL_evd: 146.4023
## LL_evd: 146.4841
## LL_evd: 146.2723
## LL_evd: 145.7731
## LL_evd: 146.5053
## LL_evd: 147.0276
## LL_evd: 146.1182
## LL_evd: 146.1479
## LL_evd: 146.5012
## LL_evd: 146.8602
## LL_evd: 146.0667
## LL_evd: 146.4012
## LL_evd: 146.2801
## LL_evd: 146.1028
## LL_evd: 146.3327
## LL_evd: 146.0966
## LL_evd: 146.4254
## LL_evd: 145.901
## LL_evd: 146.0273
## LL_evd: 146.5183
## LL_evd: 145.9801
## LL_evd: 146.8058
## LL_evd: 146.9437
## LL_evd: 146.1553
## LL_evd: 146.312
## LL_evd: 146.4264
## LL_evd: 146.9354
## LL_evd: 146.0721
## LL_evd: 145.9243
## LL_evd: 146.0956
## LL_evd: 146.6961
## LL_evd: 146.3543
## LL_evd: 146.7713
## LL_evd: 145.9033
## LL_evd: 146.0924
## LL_evd: 146.1255
## LL_evd: 145.6633
## LL_evd: 145.7914
## LL_evd: 146.1834
## LL_evd: 146.2967
## LL_evd: 146.8684
## LL_evd: 145.9435
## LL_evd: 147.0378
## LL_evd: 145.9505
## LL_evd: 146.5652
## LL_evd: 146.4105
## LL_evd: 146.228
## LL_evd: 146.7556
## LL_evd: 146.5536
## LL_evd: 146.1621
## LL_evd: 146.2329
## LL_evd: 145.6885
## LL_evd: 145.9952
## LL_evd: 146.4323
## LL_evd: 146.3742
## LL_evd: 146.0876
## LL_evd: 146.1592
## LL_evd: 145.9662
## LL_evd: 146.8657
## LL_evd: 146.0601
## LL_evd: 146.5336
## LL_evd: 146.5757
## LL_evd: 146.3253
## LL_evd: 146.4013
## LL_evd: 145.9997
## LL_evd: 146.5871
## LL_evd: 145.7387
## LL_evd: 146.6678
## LL_evd: 146.3975
## LL_evd: 145.5918
## LL_evd: 146.5317
## LL_evd: 146.9231
## LL_evd: 146.202
## LL_evd: 146.5967
## LL_evd: 146.6206
## LL_evd: 145.98
## LL_evd: 146.3342
## LL_evd: 146.2368
## LL_evd: 146.0023
## LL_evd: 146.6733
## LL_evd: 146.7509
## LL_evd: 146.0851
## LL_evd: 146.6349
## LL_evd: 146.2237
## LL_evd: 146.3211
## LL_evd: 146.2473
## LL_evd: 146.2085
## LL_evd: 146.3898
## LL_evd: 145.6857
## LL_evd: 146.135
## LL_evd: 146.4817
## LL_evd: 146.7336
## LL_evd: 146.6898
## LL_evd: 146.3134
## LL_evd: 146.6731
## LL_evd: 145.9373
## LL_evd: 146.211
## LL_evd: 146.5798
## LL_evd: 145.7075
## LL_evd: 146.9538
## LL_evd: 146.6835
## LL_evd: 146.3825
## LL_evd: 146.1211
## LL_evd: 145.8833
## LL_evd: 146.7156
## LL_evd: 146.0071
## LL_evd: 146.6806
## LL_evd: 146.2426
## LL_evd: 146.2227
## LL_evd: 146.453
## LL_evd: 146.1258
## LL_evd: 145.8465
## LL_evd: 146.6924
## LL_evd: 146.3021
## LL_evd: 147.0157
## LL_evd: 146.6116
## LL_evd: 146.555
## LL_evd: 146.3924
## LL_evd: 146.6142
## LL_evd: 146.4166
## LL_evd: 146.8208
## LL_evd: 146.0052
## LL_evd: 146.4332
## LL_evd: 145.6676
## LL_evd: 146.2173
## LL_evd: 146.0538
## LL_evd: 146.686
## LL_evd: 146.7289
## LL_evd: 145.9858
## LL_evd: 146.6508
## LL_evd: 146.4728
## LL_evd: 146.6481
## LL_evd: 146.5547
## LL_evd: 146.3759
## LL_evd: 145.766
## LL_evd: 146.7672
## LL_evd: 146.5503
## LL_evd: 146.1542
## LL_evd: 146.244
## LL_evd: 146.2482
## LL_evd: 145.7803
## LL_evd: 146.0842
## LL_evd: 146.4304
## LL_evd: 146.4013
## LL_evd: 147.0703
## LL_evd: 146.7642
## LL_evd: 146.4276
## LL_evd: 146.372
## LL_evd: 146.2377
## LL_evd: 146.6976
## LL_evd: 146.3009
## LL_evd: 146.0937
## LL_evd: 146.4618
## LL_evd: 145.9014
## LL_evd: 146.8597
## LL_evd: 145.9604
## LL_evd: 146.7366
## LL_evd: 146.8254
## LL_evd: 146.4599
## LL_evd: 146.4137
## LL_evd: 146.2719
## LL_evd: 146.6703
## LL_evd: 146.2802
## LL_evd: 146.3497
## LL_evd: 146.5821
## LL_evd: 146.5752
## LL_evd: 147.112
## LL_evd: 146.6952
## LL_evd: 146.1569
## LL_evd: 146.2965
## LL_evd: 146.0582
## LL_evd: 146.2527
## LL_evd: 146.3712
## LL_evd: 146.3832
## LL_evd: 146.3701
## LL_evd: 146.1412
## LL_evd: 146.4243
## LL_evd: 146.7697
## LL_evd: 146.5345
## LL_evd: 146.8702
## LL_evd: 146.5254
## LL_evd: 146.1377
## LL_evd: 146.8953
## LL_evd: 146.5868
## LL_evd: 146.5516
## LL_evd: 146.8515
## LL_evd: 146.0712
## LL_evd: 146.3356
## LL_evd: 146.0733
## LL_evd: 145.8681
## LL_evd: 146.8401
## LL_evd: 146.3198
## LL_evd: 145.9002
## LL_evd: 146.7556
## LL_evd: 146.2747
## LL_evd: 146.3172
## LL_evd: 146.1007
## LL_evd: 146.731
## LL_evd: 146.2238
## LL_evd: 146.4706
## LL_evd: 146.0216
## LL_evd: 146.6817
## LL_evd: 145.9051
## LL_evd: 146.2644
## LL_evd: 146.3267
## LL_evd: 146.2777
## LL_evd: 146.706
## LL_evd: 146.7326
## LL_evd: 146.74
## LL_evd: 146.1363
## LL_evd: 146.295
## LL_evd: 146.6483
## LL_evd: 146.2876
## LL_evd: 146.3252
## LL_evd: 146.7284
## LL_evd: 146.706
## LL_evd: 146.7284
## LL_evd: 146.1308
## LL_evd: 147.0856
## LL_evd: 146.6981
## LL_evd: 146.1708
## LL_evd: 146.3729
## LL_evd: 146.5145
## LL_evd: 145.8455
## LL_evd: 146.7767
## LL_evd: 146.866
## LL_evd: 146.5192
## LL_evd: 145.7926
## LL_evd: 146.0985
## LL_evd: 146.1276
## LL_evd: 146.2481
## LL_evd: 145.7679
## LL_evd: 145.7646
## LL_evd: 147.1123
## LL_evd: 146.3457
## LL_evd: 146.4711
## LL_evd: 146.1314
## LL_evd: 145.8013
## LL_evd: 146.3894
## LL_evd: 146.1882
## LL_evd: 146.4711
## LL_evd: 147.1278
## LL_evd: 145.9373
## LL_evd: 146.6322
## LL_evd: 145.6542
## LL_evd: 146.3121
## LL_evd: 146.1495
## LL_evd: 146.6394
## LL_evd: 145.9554
## LL_evd: 146.9039
## LL_evd: 146.5829
## LL_evd: 146.577
## LL_evd: 146.873
## LL_evd: 146.087
## LL_evd: 146.7411
## LL_evd: 146.8009
## LL_evd: 146.3135
## LL_evd: 146.3522
## LL_evd: 146.6487
## LL_evd: 146.577
## LL_evd: 146.6487
## LL_evd: 146.4541
## LL_evd: 146.5775
## LL_evd: 146.1567
## LL_evd: 146.2611
## LL_evd: 146.7709
## LL_evd: 146.7657
## LL_evd: 146.911
## LL_evd: 146.0957
## LL_evd: 146.5004
## LL_evd: 146.4313
## LL_evd: 146.4278
## LL_evd: 146.0447
## LL_evd: 146.3296
## LL_evd: 146.5321
## LL_evd: 146.4024
## LL_evd: 146.7316
## LL_evd: 146.4313
## LL_evd: 146.7316
## LL_evd: 146.7742
## LL_evd: 146.4124
## LL_evd: 146.3108
## LL_evd: 146.1469
## LL_evd: 146.387
## LL_evd: 146.6494
## LL_evd: 146.4355
## LL_evd: 146.011
## LL_evd: 146.0111
## LL_evd: 146.1954
## LL_evd: 146.3508
## LL_evd: 146.2439
## LL_evd: 145.9742
## LL_evd: 146.1533
## LL_evd: 146.5554
## LL_evd: 146.0495
## LL_evd: 146.69
## LL_evd: 146.043
## LL_evd: 145.9176
## LL_evd: 146.5477
## LL_evd: 146.836
## LL_evd: 146.7846
## LL_evd: 146.2173
## LL_evd: 146.578
## LL_evd: 146.4481
## LL_evd: 146.7995
## LL_evd: 146.8473
## LL_evd: 146.5903
## LL_evd: 146.1601
## LL_evd: 146.5445
## LL_evd: 146.3542
## LL_evd: 146.5442
## LL_evd: 145.7533
## LL_evd: 146.7746
## LL_evd: 146.6767
## LL_evd: 146.0184
## LL_evd: 145.8248
## LL_evd: 146.7502
## LL_evd: 146.8152
## LL_evd: 145.6973
## LL_evd: 146.314
## LL_evd: 146.2191
## LL_evd: 146.1943
## LL_evd: 146.1384
## LL_evd: 146.1217
## LL_evd: 146.4067
## LL_evd: 146.6289
## LL_evd: 147.034
## LL_evd: 146.2245
## LL_evd: 146.2957
## LL_evd: 146.9319
## LL_evd: 146.7019
## LL_evd: 146.4258
## LL_evd: 146.3568
## LL_evd: 146.5138
## LL_evd: 146.6101
## LL_evd: 146.7019
## LL_evd: 146.6101
## LL_evd: 145.9835
## LL_evd: 146.184
## LL_evd: 146.2015
## LL_evd: 146.2281
## LL_evd: 146.4868
## LL_evd: 146.6176
## LL_evd: 147.0724
## LL_evd: 146.876
## LL_evd: 146.9879
## LL_evd: 146.5324
## LL_evd: 146.3994
## LL_evd: 146.1332
## LL_evd: 146.4456
## LL_evd: 146.2969
## LL_evd: 146.2028
## LL_evd: 146.634
## LL_evd: 146.1316
## LL_evd: 146.2374
## LL_evd: 146.1895
## LL_evd: 146.4844
## LL_evd: 146.6841
## LL_evd: 146.3113
## LL_evd: 146.6946
## LL_evd: 146.7448
## LL_evd: 145.7637
## LL_evd: 147.0488
## LL_evd: 146.223
## LL_evd: 146.5651
## LL_evd: 146.4716
## LL_evd: 146.3672
## LL_evd: 146.2947
## LL_evd: 146.4372
## LL_evd: 146.646
## LL_evd: 146.9845
## LL_evd: 146.4372
## LL_evd: 146.4528
## LL_evd: 146.4676
## LL_evd: 146.2623
## LL_evd: 146.571
## LL_evd: 146.313
## LL_evd: 146.6407
## LL_evd: 145.8401
## LL_evd: 146.1942
## LL_evd: 146.9439
## LL_evd: 146.6207
## LL_evd: 146.6388
## LL_evd: 146.2989
## LL_evd: 145.8401
## LL_evd: 147.0347
## LL_evd: 146.1534
## LL_evd: 146.4038
## LL_evd: 146.8729
## LL_evd: 147.0347
## LL_evd: 146.8729
## LL_evd: 146.6388
## LL_evd: 145.8401
## LL_evd: 146.4242
## LL_evd: 146.4757
## LL_evd: 146.2989
## LL_evd: 146.2989
## LL_evd: 145.9108
## LL_evd: 145.8401
## LL_evd: 145.8957
## LL_evd: 146.2771
## LL_evd: 146.4528
## LL_evd: 146.6388
## LL_evd: 146.6388
## LL_evd: 145.8401
## LL_evd: 145.9108
## LL_evd: 145.8401
## LL_evd: 145.9108
## LL_evd: 146.2771
## LL_evd: 146.6388
## LL_evd: 146.2771
## LL_evd: 146.2771
## LL_evd: 145.8401
## LL_evd: 146.2771
## LL_evd: 145.8401
## LL_evd: 146.3979
## LL_evd: 145.8401
## LL_evd: 146.3979
## LL_evd: 145.8401
## LL_evd: 145.8401
## LL_evd: 145.8401
## LL_evd: 145.8401
## LL_evd: 145.8401
## LL_evd: 145.8401
## LL_evd: 145.8401
## LL_evd: 145.8401
## LL_evd: 145.8401
## LL_evd: 145.5339
## LL_evd: 145.5339
## LL_evd: 152.0155
## LL_evd: 156.689
## LL_evd: 146.586
## LL_evd: 148.8828
## LL_evd: 145.5121
## LL_evd: 146.6555
## LL_evd: 145.3858
## LL_evd: 145.9569
## LL_evd: 150.5588
## LL_evd: 160.6603
## LL_evd: 145.8275
## LL_evd: 150.2356
## LL_evd: 145.0946
## LL_evd: 147.2194
## LL_evd: 145.1633
## LL_evd: 146.2158
## LL_evd: 147.9703
## LL_evd: 148.4921
## LL_evd: 147.037
## LL_evd: 147.2607
## LL_evd: 146.6123
## LL_evd: 146.7517
## LL_evd: 146.7591
## LL_evd: 146.2421
## LL_evd: 147.3167
## LL_evd: 149.222
## LL_evd: 146.4098
## LL_evd: 147.1139
## LL_evd: 146.4805
## LL_evd: 146.7964
## LL_evd: 146.1111
## LL_evd: 146.9822
## LL_evd: 172.9559
## LL_evd: 192.7039
## LL_evd: 151.0708
## LL_evd: 158.3854
## LL_evd: 146.1524
## LL_evd: 149.4964
## LL_evd: 145.2895
## LL_evd: 146.9225
## LL_evd: 154.563
## LL_evd: 156.1726
## LL_evd: 147.8038
## LL_evd: 149.9558
## LL_evd: 146.526
## LL_evd: 147.7021
## LL_evd: 146.5888
## LL_evd: 146.601
## LL_evd: 152.5633
## LL_evd: 159.2938
## LL_evd: 146.9223
## LL_evd: 151.1351
## LL_evd: 146.0818
## LL_evd: 148.2222
## LL_evd: 146.3606
## LL_evd: 146.8472
## LL_evd: 153.5832
## LL_evd: 158.0237
## LL_evd: 147.2321
## LL_evd: 149.9898
## LL_evd: 146.3672
## LL_evd: 147.8138
## LL_evd: 145.9496
## LL_evd: 147.4059
## LL_evd: 151.6417
## LL_evd: 161.2935
## LL_evd: 146.363
## LL_evd: 151.2068
## LL_evd: 145.9165
## LL_evd: 148.3381
## LL_evd: 145.7235
## LL_evd: 147.6587
## LL_evd: 150.7288
## LL_evd: 151.7192
## LL_evd: 146.6409
## LL_evd: 147.9792
## LL_evd: 146.3301
## LL_evd: 146.7321
## LL_evd: 145.0179
## LL_evd: 146.9116

proc.time()-ptm

##    user  system elapsed 
## 4801.57   19.70 4825.88

Display the estimated coefficients and the standard errors. Remember the true coefficients were b0=-1.5, b1=3, rho=0.5, phi=0.25.

coef<-res$par
coef

## [1] -1.2073657  2.4153792  0.4301779  0.2427439

ses<-sqrt(diag(solve(res$hessian)))
ses

## [1] 0.095245560 0.179362321 0.005741511 0.003145638

Estimate Spatial Probit

Olga Chyzh

May 8, 2018