ちょっと特異なデータを、そこそこ正しく推定するための検討用です。結論は、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)