Y<-1:5
Y
mu<-matrix(1:100,100,100)/100*4+1
sigma2<-t(matrix(1:100,100,100)/100*6)
L<-matrix(1,100,100)
for(i in 1:5) L<-L*exp(-(Y[i]-mu)^2/sigma2/2)/sqrt(2*pi*sigma2)
contour(mu[,1],sigma2[1,],L)

##############

n<-117
m<-t(matrix(scan("c:/keith/teaching/678/fmri.txt"),3,n))
Y<-m[,1]
hot<-m[,2]
warm<-m[,3]
time<-(0:(n-1))*3

trend<-((time-mean(time))/max(time)*2)
drift<-cbind(trend^0, trend^1, trend^2, trend^3)
par(mfrow=c(2,2))
matplot(time,cbind(Y-34,lowess(time,Y)$y-34,drift,hot,warm),type="l",lwd=2)

summary(lm(Y ~ drift -1 + hot + warm))
hpw<-(hot+warm)/2
hmw<-(hot-warm)/2
summary(lm(Y ~ drift -1 + hpw + hmw))

plot(acf(lm(Y ~ drift -1 + hpw + hmw)$residuals),lwd=2)

X<-cbind(drift,hpw,hmw)

i<-matrix(1:n,n,n)
lag<-abs(i-t(i))
lag[1:5,1:5]

theta<-c(1,0)
k<-length(theta)
dV<-array(0,c(n,n,k))
E<-matrix(0,k,1)
M<-matrix(0,k,k)
dV<-array(0,c(n,n,k))
p<-dim(X)[2]
p

for(iter in 1:5)
{
V <-theta[1]*theta[2]^lag
dV[,,1]<-theta[2]^lag
dV[,,2]<-theta[1]*lag*(theta[2]+(lag==0))^(lag-1)
Vinv<-solve(V)
varbeta<-solve(t(X) %*% Vinv %*% X)
G<-varbeta %*% t(X) %*% Vinv
R<-Vinv-Vinv %*% X %*% G
m2loglr<-t(Y) %*% R %*% Y + log(det(V)) - log(det(varbeta)) - log(det(t(X) %*% X)) + (n-p)*log(2*pi)
print(m2loglr)
for(j in 1:k) E[j,1]<-sum(diag(R %*% dV[,,j])) - t(Y) %*% R %*% dV[,,j] %*% R %*% Y
print(E)
for(j1 in 1:k) for(j2 in 1:k) M[j1,j2]<-sum(diag(R %*% dV[,,j1] %*% R %*% dV[,,j2]))
Vtheta<-2*solve(M)
theta<-theta-Vtheta %*% E /2
print(theta)
}

round(V[1:5,1:5],3)
Vmh<-t(solve(chol(V)))
round(Vmh[1:5,1:5],3)
Ystar<-Vmh %*% Y
driftstar<-Vmh %*% drift
hpwstar<-Vmh %*% hpw
hmwstar<-Vmh %*% hmw
summary(lm(Ystar ~ driftstar -1 + hpwstar + hmwstar))

Yhat<-X %*% lm(Ystar ~ driftstar -1 + hpwstar + hmwstar)$coefficients
matplot(time,cbind(Y-34,lowess(time,Y)$y-34,drift,hot,warm,Yhat-34),type="l",lwd=2)

Vtheta

dvaref<-matrix(0,k,p)
for(j in 1:k) dvaref[j,]<-diag(G %*% dV[,,j] %*% t(G))
varvaref<-diag(t(dvaref) %*% Vtheta %*% dvaref)
varef<-diag(varbeta)
df<-2*varef^2/varvaref
df

##############

Vinv<-solve(V)
varbeta<-solve(t(X) %*% Vinv %*% X)
round(varbeta,5)
contrast<-c(0,0,0,0,0,1)
varef<-t(contrast) %*% varbeta %*% contrast
sqrt(varef)
cG<-t(contrast) %*% varbeta %*% t(X) %*% Vinv
ef<-cG %*% Y
ef
T<-ef/sqrt(varef)
T
dvaref<-matrix(0,k,1)
for(j in 1:k) dvaref[j,1]<-cG %*% dV[,,j] %*% t(cG)
varvaref<-t(dvaref) %*% Vtheta %*% dvaref
df<-2*varef^2/varvaref
df