extract AIC from a fixed effect spatial panel model estimation result

Question

I am using the package 'splm' of Millo and Piras(2012) to estimate a SPatial Durbin model with country (individual) fixed effect.

 ks = log(spa.sak$index_ret+2) ~ log(spa.sak$lagindex+2) + log(spa.sak$cred_def_rate+2) +   log(spa.sak$ch_in_ex_rate+2) +  
            log(spa.sak$un_inf+2) + log(spa.sak$gdp_growth+2) + log(spa.sak$cha_in_int_rate+2) + country

    ERVM<- spml(formula = ks, data = spa.sak, index = c('id','months'), listw = ERV, model = "within", lag = TRUE, effect = "individual", spatial.error = "none") 
> summary(ERVM)
Spatial panel fixed effects lag model


Call:
spml(formula = ks, data = spa.sak, index = c("id", "months"), 
    listw = ERV, model = "within", effect = "individual", lag = TRUE, 
    spatial.error = "none")

Residuals:
     Min.   1st Qu.    Median   3rd Qu.      Max. 
-0.783343 -0.028334 -0.003056  0.020040  1.199906 

Spatial autoregressive coefficient:
         Estimate Std. Error t-value Pr(>|t|)
lambda -0.0040442  0.0061785 -0.6546   0.5127

Coefficients:
                                    Estimate  Std. Error  t-value  Pr(>|t|)    
log(spa.sak$lagindex + 2)        -2.5577e-01  1.6593e-02 -15.4142 < 2.2e-16 ***
log(spa.sak$cred_def_rate + 2)    1.3427e-02  2.3244e-02   0.5777    0.5635    
log(spa.sak$ch_in_ex_rate + 2)   -5.1946e-04  2.2338e-03  -0.2325    0.8161    
log(spa.sak$un_inf + 2)           1.1360e-01  2.2112e-01   0.5138    0.6074    
log(spa.sak$gdp_growth + 2)       3.1093e+00  2.1339e+00   1.4571    0.1451    
log(spa.sak$cha_in_int_rate + 2) -1.4148e+01  3.1507e+00  -4.4905 7.105e-06 ***
countryEGYPT                      1.0659e-02  1.7940e-02   0.5942    0.5524    
countryGHANA                      2.3214e-02  2.0567e-02   1.1287    0.2590    
countryKENYA                      1.6923e-02  1.9599e-02   0.8634    0.3879    
countryMALAWI                     1.0190e-02  3.2028e-02   0.3182    0.7504    
countryMAURITIUS                  2.5843e-03  1.3530e-02   0.1910    0.8485    
countryMOROCCO                    4.2243e-03  1.4848e-02   0.2845    0.7760    
countryNAMIBIA                   -1.0416e-02  1.4408e-02  -0.7229    0.4697    
countryNIGERIA                    1.2794e-02  1.8327e-02   0.6981    0.4851    
countrySOUTH AFRICA              -9.4671e-04  1.3665e-02  -0.0693    0.9448    
countryTANZANIA                   1.9698e-02  1.9822e-02   0.9937    0.3204    
countryTUNISIA                    1.4579e-02  1.4895e-02   0.9788    0.3277    
countryUGANDA                     3.1850e-02  2.1042e-02   1.5137    0.1301    
countryZAMBIA                     3.3671e-02  2.1202e-02   1.5881    0.1123    
countryZIMBABWE                   1.4299e-02  3.1798e-02   0.4497    0.6529    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 > str(ERVM)
    List of 15
     $ coefficients : Named num [1:21] -0.004044 -0.255767 0.013427 -0.000519 0.113601 ...
      ..- attr(*, "names")= chr [1:21] "lambda" "log(spa.sak$lagindex + 2)" "log(spa.sak$cred_def_rate + 2)" "log(spa.sak$ch_in_ex_rate + 2)" ...
     $ errcomp      : NULL
     $ vcov         : num [1:21, 1:21] 3.82e-05 0.00 0.00 0.00 0.00 ...
     $ spat.coef    : Named num -0.00404
      ..- attr(*, "names")= chr "lambda"
     $ vcov.errcomp : NULL
     $ residuals    : Named num [1:3390] -0.1195 0.0395 0.0151 -0.0186 0.0285 ...
      ..- attr(*, "names")= chr [1:3390] "5-01/01/2000" "5-01/01/2001" "5-01/01/2002" "5-01/01/2003" ...
     $ fitted.values: num [1:3390, 1] 0.681 0.712 0.704 0.725 0.719 ...
      ..- attr(*, "dimnames")=List of 2
      .. ..$ : chr [1:3390] "5-01/01/2000" "5-01/01/2001" "5-01/01/2002" "5-01/01/2003" ...
      .. ..$ : NULL
      ..- attr(*, "names")= chr [1:3390] "5-01/01/2000" "5-01/01/2001" "5-01/01/2002" "5-01/01/2003" ...
     $ sigma2       : num 0.0192
     $ type         : chr "fixed effects lag"
     $ model        :'data.frame':  3390 obs. of  21 variables:
      ..$ y                               : num [1:3390] 0.537 0.731 0.697 0.685 0.728 ...
      ..$ log.spa.sak.lagindex...2.       : num [1:3390] 0.851 0.71 0.752 0.674 0.731 ...
      ..$ log.spa.sak.cred_def_rate...2.  : num [1:3390] 2.83 2.83 2.83 2.83 2.83 ...
      ..$ log.spa.sak.ch_in_ex_rate...2.  : num [1:3390] 0.696 0.69 0.711 0.768 0.695 ...
      ..$ log.spa.sak.un_inf...2.         : num [1:3390] 0.697 0.694 0.693 0.692 0.694 ...
      ..$ log.spa.sak.gdp_growth...2.     : num [1:3390] 0.693 0.693 0.693 0.693 0.693 ...
      ..$ log.spa.sak.cha_in_int_rate...2.: num [1:3390] 0.693 0.693 0.693 0.693 0.693 ...
      ..$ countryEGYPT                    : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryGHANA                    : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryKENYA                    : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryMALAWI                   : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryMAURITIUS                : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryMOROCCO                  : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryNAMIBIA                  : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryNIGERIA                  : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countrySOUTH.AFRICA             : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryTANZANIA                 : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryTUNISIA                  : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryUGANDA                   : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryZAMBIA                   : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
      ..$ countryZIMBABWE                 : num [1:3390] 0 0 0 0 0 0 0 0 0 0 ...
     $ call         : language spml(formula = ks, data = spa.sak, index = c("id", "months"), listw = ERV, model = "within", effect = "individual| __truncated__
     $ logLik       : num 1895
     $ method       : chr "eigen"
     $ effects      : chr "spfe"
     $ res.eff      :List of 2
      ..$         :List of 7
      .. ..$ res.sfe   : num [1:15, 1] 0.00264 -0.00152 0.0017 0.0023 0.0106 ...
      .. .. ..- attr(*, "dimnames")=List of 2
      .. .. .. ..$ : chr [1:15] "1" "2" "3" "4" ...
      .. .. .. ..$ : NULL
      .. ..$ res.se.sfe: Named num [1:15] 2.55 2.55 2.55 2.55 2.55 ...
      .. .. ..- attr(*, "names")= chr [1:15] "1" "2" "3" "4" ...
      .. ..$ intercept : num [1, 1] 8.43
      .. ..$ res.se.con: num [1, 1] 2.55
      .. ..$ xhat      : num [1:3390, 1] 0.681 0.712 0.704 0.725 0.719 ...
      .. .. ..- attr(*, "dimnames")=List of 2
      .. .. .. ..$ : chr [1:3390] "5-01/01/2000" "5-01/01/2001" "5-01/01/2002" "5-01/01/2003" ...
      .. .. .. ..$ : NULL
      .. ..$ N.vars    : int 35
      .. ..$ res.e     : num [1:3390, 1] -0.1195 0.0395 0.0151 -0.0186 0.0285 ...
      .. .. ..- attr(*, "dimnames")=List of 2
      .. .. .. ..$ : chr [1:3390] "5-01/01/2000" "5-01/01/2001" "5-01/01/2002" "5-01/01/2003" ...
      .. .. .. ..$ : NULL
      ..$ res.corr: NULL
     - attr(*, "class")= chr "splm"

I am successful to obtain estimation results except the following two.

1) As it is shown above, I am only able to get values for \lambda and \beta coefficients. I can not get the \theta estimation result or value that relates with the spatially lagged explanatory variables which shows the spatio-temporal relationships. How can extract it?

2) I am not able to extract AIC and BIC for model comparison. In the first case, my trials weren't successful, and For the second case, I have tried

AICsplm(ERVMX, criterion="AIC")

But I got an error.

How can I fix these issues?

UPDATE: I have solved the issue and please refer below.

Answer 1

Construct the spatial lag for both response and predictor variable in advance and add in to the formula and run.

As far as AIC is concerned, there is a need to download AICsplm.R from github download AICsplm.R to run the code.

source(AICsplm.R) 
AICsplm(ERVM)

Answer 2

Based on the answer in this question, this function also works for splm objects:

AIC_adj <- function(mod){
  # Number of observations
  n.N   <- nrow(mod$model)
  # Residuals vector
  u.hat <- residuals(mod)
  # Variance estimation
  s.sq  <- log( (sum(u.hat^2)/(n.N)))
  # Number of parameters (incl. constant) + one additional for variance estimation
  p     <-  length(coef(mod)) + 1

  # Note: minus sign cancels in log likelihood
  aic <- 2*p  +  n.N * (  log(2*pi) + s.sq  + 1 ) 

  return(aic)
}