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 )。これをより良い方法で行う方法がわかりません。グーグルで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)

アップデート：

Vincent のコードを試してみました。私は pcse2 と vcovBK コードを試しましたが、どちらも機能しました (vcocBK から得られる相関行列をどうするかはわかりませんが)。

ただし、再分析している論文の回帰係数の推定値を再現するのにまだ問題があります。私はできる限り彼らのレシピに従っています。私が見逃している唯一のステップは、Stata で「自己相関: 一般的な AR(1)」が行われる部分だと思います。私が分析している論文には次のように書かれています。

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パッケージをロードし、関数をロードpcse2して、pcseドキュメントの指示に従って使用するだけです。pcseIMHO、以下に貼り付けられた機能は、人々によって提供されたものよりもクリーンで、柔軟性があり、堅牢です。単純なベンチマークは、私のバージョンがそれらのバージョンよりも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 に答える