モンテカルロ・シミュレーション（連続ランダムウォーク）

Question

y <- cumsum(rnorm(100,0,1)) # random normal, with small (1.0) drift. 
y.ts <- ts(y) 
x <- cumsum(rnorm(100,0,1)) 
x 
x.ts <- ts(x) 
ts.plot(y.ts,ty= "l", x.ts) # plot the two random walks 


Regression.Q1 = lm(y~x) ; summary(lm2) 
summary(Regression.Q1) 

t.test1 <- (summary(Regression.Q1)$coef[2,3]) # T-test computation 


y[t] = y[t-1] + epsilon[t] 
epsilon[t] ~ N(0,1) 
set.seed(1) 
t=1000 
epsilon=sample(c(-1,1), t, replace = 1) # Generate k random walks across time {0, 1, ... , T} 


N=T=1e3 
y=t(apply(matrix(sample(c(-1,1),N*T,rep=TRUE),ncol=T),1,cumsum)) 
y[1]<-0 
for (i in 2:t) { 
    y[i]<-y[i-1]+epsilon[i] 
} 

私がする必要が統計を保存します。 S =（t-test1、t-test2、...、t-test1000）という一連の1千t検定があります。 1000回のt検定の絶対値> 1.96、5％の有意水準での臨界値の数を数えます。シリーズがI（0）なら、あなたは約5％を見つけたでしょう。ここではそうではありません（偽回帰）。

それぞれの係数を保存するには何を追加する必要がありますか？

Answer 1

y[t] = y[t-1] + epsilon[t]に関連する投稿コードは実際の動作コードではありませんが、1000 * 2のランダムウォークをすべて保存しようとしています。実際にこれを行う必要はありません。私たちはランダムウォークの実現が何であるかではなく、t-スコアだけを気にします。

このような問題では、プロシージャを何度も複製することを目指しているため、このようなプロシージャを一度に実行する関数を最初に記述すると便利です。あなたはすでにこれのために良い作業コードを持っていました。私たちは機能（plotのようなものを不要な部分を削除）でラップする必要があります。

sim <- function() { 
    y <- cumsum(rnorm(100,0,1)) 
    x <- cumsum(rnorm(100,0,1)) 
    coef(summary(lm(y ~ x)))[2,3] 
    } 

この関数には、入力を取りません。 1回の実験ではtスコアのみを返します。

ここでは、これを1000回繰り返します。これはlmとsummary.lmの両方があるので、それは、このような単純なタスクのためにあるべきよりもはるかに遅い、しばらく時間がかかります注

S <- replicate(1000, sim()) 

（必要に応じて ?replicateを読んで）私たちは、 forループを書くことができますが、機能 replicateは簡単ですスロー。もっと速い方法が後で表示されます。

今、Sは1000個の値を持つベクトルで、「1000 t検定のシーケンス」です。 を取得するには、「時間の数1,000 t検定の絶対値が> 1.96」、私たちは、結果756は私が得るだけのものです

sum(abs(S) > 1.96) 
# [1] 756 

使用することができます。シミュレーションがランダムなので、あなたは何か違うものを得るでしょう。しかし、それはいつも予想どおりかなり大きな数字になるでしょう。

---

simの高速バージョン：

fast_sim <- function() { 
    y <- cumsum(rnorm(100,0,1)) 
    x <- cumsum(rnorm(100,0,1)) 
    y <- y - mean(y) 
    x <- x - mean(x) 
    xty <- crossprod(x,y)[1] 
    xtx <- crossprod(x)[1] 
    b <- xty/xtx 
    sigma <- sqrt(sum((y - x * b)^2)/98) 
    b * sqrt(xtx) * sigma 
    } 

この関数はsummary.lmことなく、単純な線形回帰lmせず、そしてT-スコアを計算します。

S <- replicate(1000, fast_sim()) 
sum(abs(S) > 1.96) 
# [1] 778 

別の方法はcor.testを使用することです：

fast_sim2 <- function() { 
    y <- cumsum(rnorm(100,0,1)) 
    x <- cumsum(rnorm(100,0,1)) 
    unname(cor.test(x, y)[[1]]) 
    } 

S <- replicate(1000, fast_sim()) 
sum(abs(S) > 1.96) 
# [1] 775 

のベンチマークを持ってみましょう：

system.time(replicate(1000, sim())) 
# user system elapsed 
# 1.860 0.004 1.867 

system.time(replicate(1000, fast_sim())) 
# user system elapsed 
# 0.088 0.000 0.090 

system.time(replicate(1000, fast_sim2())) 
# user system elapsed 
# 0.308 0.004 0.312 

cor.testはlm + summary.lmよりもはるかに高速ですが、手動計算はさらに高速です！