2

我目前正在学习 R。我以前没有 STATA 知识。

我想重新分析在 Stata 中完成的一项研究（具有面板校正标准误差的 xtpcse 线性回归）。我在 Stata 中找不到模型或更详细的代码，也找不到如何在 R 中重写它的任何其他提示。我为 R 安装了用于计量经济学的 plm 包。就我所知。

来自 STATA 的 .do 文件的第一行复制如下（我刚刚看到它非常不可读。这是我复制 .do 内容的 txt 文件的链接：http: //dl.dropbox.com/u /4004629/This%20was%20in%20the%20.do%20file.txt）。我不知道如何以更好的方式解决这个问题。我尝试了 google-ing STATA 和 R 比较等，但没有奏效。

我想复制的研究的所有数据都在这里：

https://umdrive.memphis.edu/rblanton/public/ISQ_data

---STATA---
Group variable:   c_code                        Number of obs      =       265
Time variable:    year                          Number of groups   =        27
Panels:           correlated (unbalanced)       Obs per group: min =         3
Autocorrelation:  common AR(1)                                 avg =  9.814815
Sigma computed by pairwise selection                           max =        14
Estimated covariances      =       378          R-squared          =    0.8604
Estimated autocorrelations =         1          Wald chi2(11)      =   8321.15
Estimated coefficients     =        15          Prob > chi2        =    0.0000

------------------------------------------------------------------------------
             |           Panel-corrected
        food |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    lag_food |   .8449038    .062589    13.50   0.000     .7222316     .967576
        ciri |   -.010843   .0222419    -0.49   0.626    -.0544364    .0327504
   human_cap |   .0398406   .0142954     2.79   0.005     .0118222    .0678591
  worker_rts |  -.1132705   .0917999    -1.23   0.217    -.2931951     .066654
    polity_4 |   .0113995    .014002     0.81   0.416    -.0160439    .0388429
 market_size |   .0322474   .0696538     0.46   0.643    -.1042716    .1687665
      income |   .0382918   .0979499     0.39   0.696    -.1536865    .2302701
 econ_growth |   .0145589   .0105009     1.39   0.166    -.0060224    .0351402
   log_trade |  -.3062828   .1039597    -2.95   0.003    -.5100401   -.1025256
  fix_dollar |  -.0351874   .1129316    -0.31   0.755    -.2565293    .1861545
    fixed_xr |  -.4941214   .2059608    -2.40   0.016     -.897797   -.0904457
    xr_fluct |   .0019044   .0106668     0.18   0.858    -.0190021    .0228109
  lab_growth |   .0396278   .0277936     1.43   0.154    -.0148466    .0941022
     english |  -.1594438   .1963916    -0.81   0.417    -.5443641    .2254766
       _cons |   .4179213   1.656229     0.25   0.801    -2.828227     3.66407
-------------+----------------------------------------------------------------
         rho |   .0819359
------------------------------------------------------------------------------

. xtpcse fab_metal lag_fab_metal ciri human_cap worker_rts polity_4 market
> income econ_growth log_trade fix_dollar fixed_xr xr_fluct lab_growth
> english, pairwise corr(ar1)

更新：

我刚刚尝试了文森特的代码。我尝试了 pcse2 和 vcovBK 代码，它们都有效（尽管我不确定如何处理来自 vcocBK 的相关矩阵）。

但是，我仍然无法重现我正在重新分析的论文中回归系数的估计值。我正在尽我所能遵循他们的食谱，我认为我唯一缺少的步骤是完成 Stata“自相关：常见 AR（1）”中的部分。我正在分析的论文说：“OLS 回归使用面板校正的标准误差（Beck/Katz '95），控制每个面板内的一阶相关性（Stata 中的 corr AR1 选项）。”

如何控制 R 中每个面板内的一阶相关性？

这是我到目前为止对数据所做的：

## run lm 
res.lm <- lm(total_FDI ~ ciri + human_cap + worker_rts + polity_4 + lag_total + market_size + income + econ_growth + log_trade + fixed_xr + fix_dollar + xr_fluct + english + lab_growth, data=D)
## run pcse
res.pcse <- pcse2(res.lm,groupN="c_code",groupT="year",pairwise=TRUE)
4

2 回答 2

4

正如 Ramnath 所提到的，pcse 包将做 Stata 所做的事情xtpcse。或者，您可以使用plm 包vcovBK()中的函数。如果您选择后一个选项，请确保您使用该选项，这是 Beck & Katz (1995) 文章所建议的以及 Stata 命令实现的内容。cluster='time'

该pcse软件包运行良好，但存在一些问题使得许多直观的用户输入无法接受，尤其是在您的数据集不平衡的情况下。您可能想尝试重写我不久前编写的函数。只需按照 pcse 文档中的说明加载pcse包、加载函数并使用它。pcse2恕我直言，下面粘贴的功能比人们提供的功能更干净、更灵活、更健壮pcse。简单的基准测试还表明，我的版本可能比他们的版本快 5 到 10 倍，这对于大型数据集可能很重要。

祝你好运！

library(Matrix)
pcse2 <- function(object, groupN, groupT, pairwise=TRUE){
  ## Extract basic model info
  groupT <- tail(as.character((match.call()$groupT)), 1)
  groupN <- tail(as.character((match.call()$groupN)), 1)
  dat <- eval(parse(text=object$call$data))

  ## Sanity checks
  if(!"lm" %in% class(object)){stop("Formula object must be of class 'lm'.")}
  if(!groupT %in% colnames(dat)){stop(paste(groupT, 'was not found in data', object$call$data))}
  if(!groupN %in% colnames(dat)){stop(paste(groupN, 'was not found in data', object$call$data))}
  if(anyDuplicated(paste(dat[,groupN], dat[,groupT]))>0){stop(paste('There are duplicate groupN-groupT observations in', object$call$data))}
  if(length(dat[is.na(dat[,groupT]),groupT])>0){stop('There are missing unit indices in the data.')}
  if(length(dat[is.na(dat[,groupN]),groupN])>0){stop('There are missing time indices in the data.')}

  ## Expand model frame to include groupT, groupN, resid columns.
  f <- as.formula(object$call$formula)
  f.expanded <- update.formula(f, paste(". ~ .", groupN, groupT, sep=" + "))
  dat.pcse <- model.frame(f.expanded, dat) 
  dat.pcse$e <- resid(object)  

  ## Extract basic model info (part II)
  N <- length(unique(dat.pcse[,groupN]))
  T <- length(unique(dat.pcse[,groupT]))
  nobs <- nrow(dat.pcse)
  is.balanced <- length(resid(object)) == N * T

  ## If balanced dataset, calculate as in Beck & Katz (1995)
  if(is.balanced){
    dat.pcse <- dat.pcse[order(dat.pcse[,groupN], dat.pcse[,groupT]),]
    X <- model.matrix(f, dat.pcse)
    E <- t(matrix(dat.pcse$e, N, T, byrow=TRUE))
    Omega <- kronecker((crossprod(E) / T), Matrix(diag(1, T)) )

  ## If unbalanced and pairwise, calculate as in Franzese (1996)
  }else if(pairwise==TRUE){
    ## Rectangularize
    rectangle <- expand.grid(unique(dat.pcse[,groupN]), unique(dat.pcse[,groupT]))
    names(rectangle) <- c(groupN, groupT)
    rectangle <- merge(rectangle, dat.pcse, all.x=TRUE)
    rectangle <- rectangle[order(rectangle[,groupN], rectangle[,groupT]),]
    valid <- ifelse(is.na(rectangle$e),0,1) 
    rectangle[is.na(rectangle)] <- 0
    X <- model.matrix(f, rectangle)
    X[valid==0,1] <- 0

    ## Calculate pcse
    E <- crossprod(t(matrix(rectangle$e, N, T, byrow=TRUE)))
    V <- crossprod(t(matrix(valid, N, T, byrow=TRUE)))
    if (length(V[V==0]) > 0){stop("Error! A CS-unit exists without any obs or without any obs in a common period with another CS-unit. You must remove that unit from the data passed to pcse().")}
    Omega <-  kronecker(E/V, Matrix(diag(1, T)))

  ## If unbalanced and casewise, caluate based on largest rectangular subset of data
  }else{ 
    ## Rectangularize
    rectangle <- expand.grid(unique(dat.pcse[,groupN]), unique(dat.pcse[,groupT]))
    names(rectangle) <- c(groupN, groupT)
    rectangle <- merge(rectangle, dat.pcse, all.x=TRUE)
    rectangle <- rectangle[order(rectangle[,groupN], rectangle[,groupT]),]
    valid <- ifelse(is.na(rectangle$e),0,1) 
    rectangle[is.na(rectangle)] <- 0
    X <- model.matrix(f, rectangle)
    X[valid==0,1] <- 0

    ## Keep only years for which we have the max number of observations
    large.panels <- by(dat.pcse, dat.pcse[,groupT], nrow) # How many valid observations per year?
    if(max(large.panels) < N){warning('There is no time period during which all units are observed. Consider using pairwise estimation.')}
    T.balanced <- names(large.panels[large.panels==max(large.panels)]) # Which years have max(valid observations)?
    T.casewise <- length(T.balanced)
    dat.balanced <- dat.pcse[dat.pcse[,groupT] %in% T.balanced,] # Extract biggest rectangular subset
    dat.balanced <- dat.balanced[order(dat.balanced[,groupN], dat.balanced[,groupT]),]
    e <- dat.balanced$e

    ## Calculate pcse as in Beck & Katz (1995)
    E <- t(matrix(dat.balanced$e, N, T.casewise, byrow=TRUE))
    Omega <- kronecker((crossprod(E) / T.casewise), Matrix(diag(1, T)))
  }

  ## Finish evaluation, clean and output
  salami <- t(X) %*% Omega %*% X
  bread <- solve(crossprod(X))
  sandwich <- bread %*% salami %*% bread
  colnames(sandwich) <- names(coef(object))
  row.names(sandwich) <- names(coef(object))
  pcse <- sqrt(diag(sandwich))
  b <- coef(object)
  tstats <- b/pcse
  df <- nobs - ncol(X)
  pval <- 2*pt(abs(tstats), df, lower.tail=FALSE)
  res <- list(vcov=sandwich, pcse=pcse, b=b, tstats=tstats, df=df, pval=pval, pairwise=pairwise, 
              nobs=nobs, nmiss=(N*T)-nobs, call=match.call())
  class(res) <- "pcse"
  return(res)
}
于 2011-04-04T18:47:15.387 回答
3

查看pcse 包，它考虑了面板纠正的标准错误。您当然必须查看 STATA 中的文档以找出所做的假设并与 pcse 进行交叉检查。

于 2011-04-04T16:19:56.940 回答