入門共分散構造分析の実際 朝野・鈴木・小島隆夫先生 講談社サイエンティフィク

#共分散構造分析のパラメータ推定 p.54

data(EuStockMarkets)
dat=data.frame(EuStockMarkets)
dat=dat[1:300,]
XX=matrix(log(as.numeric(unlist(dat))),nrow=nrow(dat))
varmat=cov(XX)

predV=function(x){
a=x[1];b=x[2];sigf1=x[3];sige1=x[4];sige2=x[5];sigx3=x[6]  
A=array(0,dim=c(4,4))
A[2,1]=b;A[3,1]=b;A[3,4]=a
A=A+diag(1,4)
V=diag(c(sigf1^2,sige1^2,sige2^2,sigx3^2))
V=A%*%V%*%t(A)  
return(V)  
}

f=function(x){
V=predV(x)
return(sum(diag(solve(V)%*%varmat))-log(det(solve(V)%*%varmat)))
}



#勾配降下法

param=rep(1,6)
ite=10^(4)
eta=10^(-4)
h=0.001

for(l in 1:ite){
for(s in 1:length(param)){
param_sub=param
param_sub[s]=param_sub[s]+h
param[s]=param[s]-eta*(f(param_sub)-f(param))/h
}  
print(f(param))  
}

#予測された分散行列
predV(param)

#実際の分散行列
varmat

#適合度のカイ二乗検定
p=length(c(diag(varmat),varmat[lower.tri(varmat)]))
q=length(param)
df=p*(p+1)/2-q
kai2=(nrow(XX)-1)*f(param)
qchisq(1-0.05,df)