t→x,yの擬似相関
下の例で計算すると0.96。
線形トレンド
観測方程式
の状態空間で、x,y(変数名の方)の残差の相関をstanで計算すると-0.12程度になった。
N.t=60 mu.x = numeric(N.t) x = numeric(N.t) mu.y = numeric(N.t) x=y=c() mu.x[1] = mu.y[1] = 10 s1 = 2 s2 = 1 for (t in 1:(N.t - 1)) { x[t] = rnorm(1, mu.x[t], s1) mu.x[t + 1] = rnorm(1, mu.x[t] + 1, s2) y[t] = rnorm(1, mu.y[t], s1) mu.y[t + 1] = rnorm(1, mu.y[t] + 1, s2) } x[N.t] = rnorm(1, mu.x[N.t - 1], s1) y[N.t] = rnorm(1, mu.y[N.t - 1], s1) scodeaa=" data{ int N; real x[N]; real y[N]; } parameters{ real sigma_x; real sigma_y; real trend_x[N]; real trend_y[N]; real sigma_trend_x; real sigma_trend_y; } model{ for(t in 3:N){ trend_x[t]~normal(trend_x[t-1]*2-trend_x[t-2],sigma_trend_x); trend_y[t]~normal(trend_y[t-1]*2-trend_y[t-2],sigma_trend_y); } x~normal(trend_x,sigma_x); y~normal(trend_y,sigma_y); } " library(rstan) fit=stan(model_code = scodeaa,data = list(N=N.t,x=x,y=y),chain=1,iter=200) la=extract(fit) tr_x=apply(la$trend_x,2,mean) tr_y=apply(la$trend_y,2,mean) cor(x-tr_x,y-tr_y) #>[1] -0.1202267