餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

打ち切りリセット付単位根過程のAR(1)での推定

ちょっと特異なデータを、そこそこ正しく推定するための検討用です。結論は、AR(1)にしておけば無問題でした。

### データ生成 ###
genSample <- function(n=250){
  ex <- rnorm(n)
  x <- rep(NA, length = n)
  x[1] <- ex[1] * 0.2
  for (j in 2:n) {
# 係数0.98で微妙に単位根ではないですが、250だとほぼ単位根過程と言うことで
    x[j] <- x[(j - 1)] * 0.98 + ex[j] * 0.2
  }
  x <- x + 5
  
  ey<-rnorm(n)
  y<-rep(NA,length=n)
  y[1]<-60+ey[1]*5
  
  d <- numeric(n)
  dc <- 0
  
  g <- numeric(n)
  gc <- 0
  
  chg <- numeric(n)
  
  for(i in 2:n){
    if(y[(i-1)]>20){
      y[i]<-y[(i-1)]-2+ey[i]*5
      dc <- dc + 1
      chg[i] <- 0
    }
    else{
      y[i]<-60+ey[i]*5
      dc <- 0
      gc <- gc + 1
      chg[i] <- 1
    }
    d[i] <- dc
    g[i] <- gc
  }
  data.frame(
    y2 = - y[-1],
    y1 = y[-length(y)],
    x = x[-1],
    d = d[-1],
    g = g[-1],
    chg = chg[-1]
  )
}

### m回施行して傾向を見る ###
n <- 250
m <- 3000
Pr_LVL <- numeric(m)
Pr_AR1 <- numeric(m)
for(c in 1:m){
  df1 <- genSample(n)

  r1 <- lm(y2 ~ x + factor(g) + d, data=df1)
  Pr_LVL[c] <- coef(summary(r1))[,4][2]

  r2 <- lm(y2 ~ y1 + x + factor(g) + chg, data=df1)
  Pr_AR1[c] <- coef(summary(r2))[,4][3]
}

hist(Pr_LVL, breaks=10)
hist(Pr_AR1, breaks=10)