Py学习  »  Python

推荐&空间计量及Python、Stata、R操作代码

计量经济学服务中心 • 3 年前 • 475 次点击  

来源:https://geodacenter.github.io/GeoDaSpace/


本文将主要介绍使用空间计量Python、Stata、R软件来复制实现

Technical Aspects of Implementing GMM Estimation of the Spatial Error Model in PySAL and GeoDaSpace∗ by Luc Anselin, Pedro V. Amaral and Daniel Arribas-Bel中空间误差模型的GMM估计的技术方面的论文结果。


1

Stata代码汇总



代码为

/* This routine assumes the following files are in the current Stata working directory: *//* - NAT.shp *//* - NAT.dbf */
/* Create row-standardized Queen weights matrix */shp2dta using NAT, database(nat) coordinates(natxy) genid(ids) gencentroids(nat_c) replaceuse natspmat contiguity w using natxy, norm(row) id(ids) replace
/* Table 2 - Spatial error model with exogenous variables and homoskedasticity (GM_Error_Hom): */spivreg HR90 RD90 UE90, el(w) id(ids)
/* Table 3 - Spatial error model with endogenous variables and homoskedasticity (GM_Endog_Error_Hom): */spivreg HR90 RD90 (UE90 = FP89), el(w) id(ids)
/* Table 4 - Combo model with homoskedasticity (GM_Combo_Hom): */ spivreg HR90 RD90 UE90, dl(w) el(w) id(ids)
/* Table 5 - Spatial error model with exogenous variables and heteroskedasticity (GM_Error_Het): */spivreg HR90 RD90 UE90, el(w) id(ids) het
/* Table 6 - Spatial error model with endogenous variables and heteroskedasticity (GM_Endog_Error_Het): */spivreg HR90 RD90 (UE90 = FP89), el(w) id(ids) het
/* Table 7 - Combo model with heteroskedasticity (GM_Combo_Het): */spivreg HR90 RD90 UE90, el(w) dl(w) id(ids) het




2

Python部分



首先导入相关库,代码为

import numpy as npimport pysal as ps


然后打开数据库并创建变量作为numpy数组,代码为:

db = ps.open('nat.dbf','r')

其中,为dbf格式数据,r表示打开方式为可读




Regression 1: use all defaults, i.e., A1='het' changed by LA to A1='hom_sc' as in Drukker

相关代码为:

reg1 = ps.spreg.GM_Error_Hom(y,x,w,name_y='HR90',name_x = x_names, name_w = 'nat_queen.gal',name_ds='NAT')

然后打印上述结果

相关代码为:

print reg1.summary

结果为:

Regression 1a: use A1 matrix as in Anselin with A1='het'

相关代码为:

reg1a = ps.spreg.GM_Error_Hom(y,x,w,A1='het',name_y='HR90',name_x = x_names, name_w = 'nat_queen.gal',name_ds='NAT')print reg1a.summary

结果为:


Regression 1b: use A1 matrix as in Anselin, without scaling factor, A1='hom'

相关代码为:

reg1b = ps.spreg.GM_Error_Hom(y,x,w,A1='hom',name_y='HR90',name_x = x_names, name_w = 'nat_queen.gal',name_ds='NAT')
print reg1b.summary

结果为:

Regression 1c: use A1 with scaling factor, but force estimation as 2SLS (as in Stata)


Need to use base class to force passing a constant term, since the user class puts constant in x

相关代码为:

ones = np.ones(y.shape)reg1c = ps.spreg.error_sp_hom.BaseGM_Endog_Error_Hom(y,ones,yend=x, q=x, w=w, A1='hom_sc')print reg1c.betasprint map(np.sqrt, reg1c.vm.diagonal())

结果为:



1.2  Endogenous


建立新的外生、内生和工具变量


相关代码为:

xex = x[:,0]xex.shape = (len(hr90),1)

yend = x[:,1]yend.shape = (len(hr90),1)

fp89 = db.by_col("FP89")q = array(fp89)q.shape = (len(fp89),1)
#Error Hom model with endogeneity - defaultreg2 = ps.spreg.GM_Endog_Error_Hom(y,xex,yend,q,w,name_y="HR90",name_x=["RD90"],name_yend=["UE90"],name_q=["FP89"],\ name_w="natqueen.gal",name_ds="nat.shp")print reg2.summary

Error hom model with endogeneity -- A1 = 'het'reg2a = ps.spreg.GM_Endog_Error_Hom(y,xex,yend,q,w,A1='het',name_y="HR90",name_x=["RD90"],name_yend=["UE90"],name_q=["FP89"],\ name_w="natqueen.gal" ,name_ds="nat.shp")print reg2a.summary

结果为:


1.3  Combo Model¶


Error hom model -- combo -- instrument lag order set to 2

相关代码为:

reg3 = ps.spreg.GM_Combo_Hom(y,x,w=w,w_lags=2,name_y='HR90',name_x = x_names, name_w = 'nat_queen.gal',name_ds='NAT')
print reg3.summary


# Error hom model -- combo -- instrument lag order set to 1reg3a = ps.spreg.GM_Combo_Hom(y,x,w=w,name_y='HR90',name_x = x_names, name_w = 'nat_queen.gal',name_ds='NAT')print reg3a.summary

# Error Hom model -- combo -- instrument lag order two and 'het' optionreg3b = ps.spreg.GM_Combo_Hom(y,x,w=w,w_lags=2,A1='het',name_y='HR90',name_x = x_names, name_w = 'nat_queen.gal',name_ds='NAT')print reg3b.summary

结果为:

2  Heteroskedasticity¶

2.1  Exogenous only¶


相关代码为:

# Defaults, two step as in Drukker et al.reg4 = ps.spreg.GM_Error_Het(y,x,w,name_y="HR90",name_x=x_names,name_w="natqueen.gal",name_ds="nat.shp")print reg4.summary
Step1c to True as in Arraiz et al (same as sphet1)reg4a = ps.spreg.GM_Error_Het(y,x,w,step1c=True,name_y="HR90" ,name_x=x_names,name_w="natqueen.gal",name_ds="nat.shp")print reg4a.summary

# Setup to mimic Stata "OLS" estimationreg4b = ps.spreg.error_sp_het.BaseGM_Endog_Error_Het(y,ones,yend=x, q=x, w=w)print reg4b.betas



结果为:

2.2  Endogenous


相关代码为:

# Standard case -- defaultsreg5 = ps.spreg.GM_Endog_Error_Het(y,xex,yend,q,w,name_y="HR90",name_x=["RD90"],name_yend=["UE90"],name_q=["FP89"],\       name_w="natqueen.gal",name_ds="nat.shp")
print reg5.summary

# step1c = Truereg5a = ps.spreg.GM_Endog_Error_Het(y,xex,yend,q,w,step1c=True,name_y="HR90",name_x=["RD90"],name_yend=["UE90"],name_q=["FP89"],\ name_w="natqueen.gal",name_ds="nat.shp")print reg5a.summary

结果为:


2.3  Combo model


相关代码为:

#  Error het  -- combo -- instrument lag order set to 2reg6 = ps.spreg.GM_Combo_Het(y,x,w=w,w_lags=2,name_y='HR90',name_x = x_names, name_w = 'nat_queen.gal',name_ds='NAT')print reg6.summary#  Error het model -- combo -- instrument lag order set to 1reg6a = ps.spreg.GM_Combo_Het(y,x,w=w,name_y='HR90',name_x = x_names, name_w = 'nat_queen.gal',name_ds='NAT')print reg6a.summary
# Error Het model -- combo -- instrument lag order two and step1c=True option
reg6c = ps.spreg.GM_Combo_Het(y,x,w=w,w_lags=2,step1c=True,name_y='HR90',name_x = x_names, name_w = 'nat_queen.gal',name_ds='NAT')print reg6c.summary

结果为:


3

R代码汇总


  1. %load_ext rmagic

GMM Comparison

This document reproduces the R results in the paper *Technical Aspects of Implementing GMM Estimation of the Spatial Error Model in PySAL and GeoDaSpace∗ by Luc Anselin, Pedro V. Amaral and Daniel Arribas-Bel.

The main library used is sphet, in its stable version (available from CRAN) and its bleeding edge development version (available from R-forge here). In this document, we use the development version as available from the master trunk on Novemer 28th., 2012.

In addition, the following dependencies are required:

  • spdep

  • sem

  • foreign

Load the dependency libraries (note loading them will produce several output lines, keep scrolling down).

  1. %%R


  2. library(spdep, quiet=T)

  3. library(sem)

  4. library(foreign)

  1. Loading required package : lattice


  2. Attaching package: lattice


  3. The following object(s) are masked from package:boot’:


  4. melanoma


  5. Loading required package: foreign

  6. Checking rgeos availability: FALSE

  7. Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,

  8. which has a restricted licence. It is disabled by default;

  9. to enable gpclib, type gpclibPermit()

  10. deldir 0.0-21

  11. Loading required package: matrixcalc

Data

To replicate the results, you will need the  NAT dataset in your machine. Throughout the document, we assume you have both the dbf file and a queen contiguity gal file created from the shapefile in your working directory. If this is the case, the path for the files should just be their names:

  1. %%R

  2. gal.path 'NAT_queen.gal'

  3. dbf.path 'NAT.dbf'

NAT_queen.gal is not included but you can easily create it in PySAL or in OpenGeoDa. Once they are both in place, load them into the session:

  1. %%R

  2. w read.gal(gal.path, override.id=TRUE)

  3. w nb2listw(w)

  4. dbf read.dbf(dbf.path)

sphet2 code

The code in this section reproduces all the results in the paper reported to come from the development version of sphet (codenamed sphet2 througout the paper).

Next, source sphet2. First specify the local path in your machine where the folder with the code is:

  1. %%R

  2. sphet2.path "/home/dani/code/R/sphet/pkg/R/"

If the files are in the same directory as this document, simply replace the path by "" (that is: sphet2.path""). Now you can load sphet2 code:

  1. %%R

  2. source(paste(sphet2.path, "circular.R", sep=''))

  3. source(paste(sphet2.path, "dist_functions.R", sep=''))

  4. source(paste(sphet2.path, "distance.R", sep=''))

  5. source(paste(sphet2.path, "gm_functions.R", sep=''))

  6. source(paste(sphet2.path, "gs2slshet.R", sep=''))

  7. source(paste(sphet2.path, "gwt2dist.R", sep=''))

  8. source(paste(sphet2.path, "ivreg.R", sep=''))

  9. source(paste(sphet2 .path, "kernelsfun.R", sep=''))

  10. source(paste(sphet2.path, "listw2dgCMatrix.R", sep=''))

  11. source(paste(sphet2.path, "Omega.R", sep=''))

  12. source(paste(sphet2.path, "s2slshac.R", sep =''))

  13. source(paste(sphet2.path, "spreg.R", sep=''))

  14. source(paste(sphet2.path, "summary.sphet.R", sep=''))

  15. source(paste(sphet2.path, "twostagels.R", sep=''))

  16. source(paste(sphet2.path, "utilities.R", sep=''))

Models

Table 2 (GMErrorHom)

  1. %%R

  2. model spreg(HR90 ~ RD90 + UE90, data = dbf,listw = w, model = "error" )

  3. print(summary( model))


  1. Generalized stsls


  2. Call:

  3. spreg(formula = HR90 ~ RD90 + UE90, data = dbf, listw = w, model = "error")


  4. Residuals:

  5. Min. 1st Qu. Median Mean 3rd Qu. Max.

  6. -18.400 -3.420 -0.662 0.019 2.530 67.100


  7. Coefficients:

  8. Estimate Std. Error t-value Pr(>|t|)

  9. (Intercept) 6.676229 0.349809 19.0853 <2e-16 ***

  10. RD90 3.944978 0.155334 25.3967 <2e-16 ***

  11. UE90 -0.077044 0.047125 -1.6349 0.1021

  12. rho 0.414905 0.019389 21.3987 <2e-16 ***

  13. ---

  14. Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 1


Table 3 (GMEndogError_Hom)

  1. %%R

  2. model spreg(HR90 ~ RD90, endog = dbf$UE90 , instruments = cbind(dbf$FP89), data = dbf,listw = w, model = "error",verbose=TRUE)

  3. print(summary(model))

  1. function: 169.7589 lambda: 0.2

  2. function: 169.7588 lambda: 0.2

  3. function: 81.68172 lambda: 0.9

  4. function: 81.67946 lambda: 0.8999915

  5. function: 47.50417 lambda: 0.7771504

  6. function: 47.50421 lambda: 0.7771506

  7. function: 1411.72 lambda: -0.2228496

  8. function: 21.65544 lambda: 0.6771504

  9. function: 21.6555 lambda: 0.6771507

  10. function: 3.706638 lambda: 0.4771504

  11. function: 3.706629 lambda: 0.4771505

  12. function: 1.78937 lambda: 0.5402429

  13. function: 1.789369 lambda: 0.5402429

  14. function: 1.441354 lambda: 0.5235599

  15. function: 1.441354 lambda: 0.5235599

  16. function: 1.438217 lambda: 0.5217685

  17. function: 1.438217 lambda: 0.5217686

  18. function: 1.438211 lambda: 0.5218404

  19. function: 1.438211 lambda: 0.5218404

  20. function: 1.438211 lambda: 0.5218402

  21. function: 1.438211 lambda: 0.5218405

  22. function: 1.438211 lambda: 0.5218398

  23. function: 1.438211 lambda: 0.5218402

  24. function: 0.004362049 lambda: 0.5218402

  25. function: 0.00436205 lambda: 0.5218402

  26. function: 0.003979579 lambda: 0.4779536

  27. function: 0.003979577 lambda: 0.4779536

  28. function: 0.00374883 lambda: 0.4952313

  29. function: 0.00374883 lambda: 0.4952313

  30. function: 0.003748192 lambda: 0.4943871

  31. function: 0.003748192 lambda: 0.4943871

  32. function: 0.003748192 lambda: 0.4943604

  33. function: 0.003748192 lambda: 0.4943605

  34. function: 0.003748192 lambda: 0.4943605

  35. function: 0.003748192 lambda: 0.4943614

  36. function: 0.003748192 lambda: 0.4943595

  37. function: 0.003748192 lambda: 0.4943605


  38. Generalized stsls


  39. Call:

  40. spreg(formula = HR90 ~ RD90, data = dbf, listw = w, endog = dbf$UE90,

  41. instruments = cbind(dbf$FP89), model = "error", verbose = TRUE)


  42. Residuals:

  43. Min. 1st Qu. Median Mean 3rd Qu. Max.

  44. -32.700 -4.100 0.141 0.034 3.950 65.400


  45. Coefficients:

  46. Estimate Std. Error t-value Pr(>|t|)

  47. (Intercept) 21.060572 1.538510 13.689 < 2.2e-16 ***

  48. RD90 8.241952 0.488835 16.860 < 2.2e-16 ***

  49. endogenous -2.243751 0.228954 -9.800 < 2.2e-16 ***

  50. rho 0.494360 0.021735 22.745 < 2.2e-16 ***

  51. ---

  52. Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 1


Table 4 (GMComboHom)

  1. %%R

  2. model spreg(HR90 ~ RD90 + UE90, data = dbf,listw = w, het = TRUE, model = "error" )

  3. print(summary(model))


  1. Generalized stsls


  2. Call:

  3. spreg(formula = HR90 ~ RD90 + UE90, data = dbf, listw = w, model = "error",

  4. het = TRUE)


  5. Residuals:

  6. Min. 1st Qu. Median Mean 3rd Qu. Max.

  7. -18.400 -3.410 -0.664 0.019 2.530 67.100


  8. Coefficients:

  9. Estimate Std. Error t-value Pr(>|t|)

  10. (Intercept) 6.658555 0.474454 14.0341 <2e-16 ***

  11. RD90 3.941664 0.259899 15.1661 <2e-16 ***

  12. UE90 -0.074496 0.061092 -1.2194 0.2227

  13. rho 0.474047 0.023665 20.0314 <2e-16 ***

  14. ---

  15. Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 1


Table 6 (GMEndogError_Het)

  1. %%R

  2. model spreg(HR90 ~ RD90, endog = dbf$UE90, instruments = dbf$FP89, data = dbf,listw = w, het = TRUE, model = "error")

  3. print (summary(model))


  1. Generalized stsls


  2. Call:

  3. spreg(formula = HR90 ~ RD90, data = dbf, listw = w, endog = dbf$UE90,

  4. instruments = dbf$FP89, model = "error", het = TRUE)


  5. Residuals:

  6. Min. 1st Qu. Median Mean 3rd Qu. Max.

  7. -32.700 -4.100 0.146 0.035 3.950 65.400


  8. Coefficients:

  9. Estimate Std. Error t-value Pr(>|t|)

  10. (Intercept) 21.028779 2.562926 8.2050 2.306e-16 ***

  11. RD90 8.237610 0.781683 10.5383 < 2.2e-16 ***

  12. endogenous -2.239152 0.390154 -5.7392 9.515e-09 ***

  13. rho 0.466672 0.029767 15.6777 < 2.2e-16 ***

  14. ---

  15. Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 1


Table 7 (GMComboHet)

  1. %%R

  2. model spreg(HR90 ~ RD90 + UE90, data = dbf,listw = w, het = TRUE, model = "sarar" )

  3. print(summary(model))


  1. Generalized stsls


  2. Call :

  3. spreg(formula = HR90 ~ RD90 + UE90, data = dbf, listw = w, model = "sarar",

  4. het = TRUE)


  5. Residuals:

  6. Min. 1st Qu. Median Mean 3rd Qu. Max.

  7. -18.500 -3.470 -0.704 0.015 2.540 67.000


  8. Coefficients:

  9. Estimate Std. Error t-value Pr(>|t|)

  10. (Intercept) 6.940607 0.859998 8.0705 7.001e-16 ***

  11. RD90 4.007396 0.326093 12.2891 < 2.2e-16 ***

  12. UE90 -0.095717 0.066417 -1.4412 0.1495

  13. lambda -0.022001 0.087638 -0.2511 0.8018

  14. rho 0.558391 0.050685 11.0168 < 2.2e-16 ***

  15. ---

  16. Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 1



Python社区是高质量的Python/Django开发社区
本文地址:http://www.python88.com/topic/62412
 
475 次点击