経時データ解析の数理　藤越 康祝先生　朝倉書店

library(dplyr)

#p.24 多群の場合

data("EuStockMarkets")

t=seq(0,1859*0.04,0.04)+rep(1991.496,1860)

t=t/1000

times=cbind(t,t^2);colnames(times)=c("t","t2")

groups=ncol(EuStockMarkets)

data=data.frame(times,EuStockMarkets)

times=as.matrix(cbind(rep(1,nrow(times)),times))

V=t(times)%*%times

sitas=svd(V)$u%*%diag(1/svd(V)$d)%*%t(svd(V)$v)%*%t(times)%*%EuStockMarkets

pre=times%*%sitas

sij2=(EuStockMarkets-pre)^2

N=prod(dim(EuStockMarkets))

sig_hat=sum(sij2)/N

sig_childer=sig_hat*N/(N-groups*ncol(times))

#回帰係数の信頼区間幅

qt(1-0.05/2,df=N-groups*ncol(times))*sqrt(diag(solve(V))*sig_childer)

#各標本の信頼区間幅

C_sita=t(sitas)%*%times[1,]

ncol(times)*qf(1-0.05,groups*ncol(times),N-groups*ncol(times))*sig_chider*t(c(times[1,]))%*%solve(V)%*%t(t(c(times[1,])))

t=c(8,8.5,9,9.5)

mat=matrix(c(47.8,46.4,46.3,45.1,47.6,52.5,51.2,49.8,48.1,45,51.2,48.5,52.1,48.2,49.6,50.7,47.2,53.3,46.2,46.3,48.8,47.3,46.8,45.3,48.5,53.2,53,50,50.8,47,51.4,49.2,52.8,48.9,50.4,51.7,47.7,54.6,47.5,47.6,49,47.7,47.8,46.1,48.9,53.3,54.3,50.3,52.3,47.3,51.6,53,53.7,49.3,51.2,52.7,48.4,55.1,48.1,51.3,49.7,48.4,48.5,47.2,49.3,53.7,54.5,52.7,54.4,48.3,51.9,55.5,55,49.8,51.8,53.3,49.5,55.3,48.4,51.8),ncol=4)

t_vec=c(rep(8,20),rep(8.5,20),rep(9,20),rep(9.5,20))

sita1=cov(c(mat),c(t_vec))/var(c(t_vec),c(t_vec))

sita0=mean(c(mat))-mean(t_vec)*sita1

pi=ncol(mat);n=nrow(mat)

s=c()

for(j in 1:n){

s=c(s,t(c(mat[j,]))%*%(diag(rep(1,pi))-t(t(t))%*%t(t)/sum(t^2))%*%t(t(c(mat[j,])))/(n-1))

}

sig_hat=sum(s)/(pi-1)