## com10a of february 2022,
#  gaussian process, parameters mu, sigma, autocorrelation rho.

#  how well does ML work? i use loglikelihood profile for rho,
#  then read off ML for mu and sigma,
#  and check in simulations how well \hat\mu, \hat\sigma, \hat\rho work.

nn <- 100
sigma0 <- 1.00
xx <- sigma0*rnorm(nn)

one <- 1 + 0*(1:nn)
dim(one) <- c(nn,1) 
A <- 0*(1:nn) %*% t(1:nn) 

Ahat <- 0*(1:nn) %*% t(1:nn) 

sim <- 100
ML <- 0*(1:sim) %*% t(1:3)

for (ss in 1:sim)
{{
xx <- sigma0*rnorm(nn)

logL <- function(rho)
{
# first the A matrix 
for (i in 1:nn)
{for (j in 1:nn)
{A[i,j] <- rho^(abs(j-i))}}
#
Ainv <- solve(A)
muhat <- t(one) %*% Ainv %*% xx / t(one) %*% Ainv %*% one
dim(muhat) <- NULL 
Q0 <- t(xx-one*muhat) %*% Ainv %*% (xx-one*muhat)
sigmahat <- sqrt(Q0/nn)
#
-nn*log(sigmahat) -0.5*log(det(A))
}

rhoval <- seq(-0.50,0.50,by=0.02)
logLval <- 0*rhoval

for (k in 1:length(rhoval))
{
logLval[k] <- logL(rhoval[k])
}

matplot(rhoval,logLval,type="l")
aux <- rhoval[logLval == max(logLval)]

minuslogL <- function(rho)
{-logL(rho)}

nils <- nlm(minuslogL,aux,hessian=T)
rhohat <- nils$estimate

for (i in 1:nn)
{for (j in 1:nn)
{Ahat[i,j] <- rhohat^(abs(j-i))}}
#
Ahatinv <- solve(Ahat)
muhat <- t(one) %*% Ahatinv %*% xx / t(one) %*% Ahatinv %*% one
dim(muhat) <- NULL 
Q0 <- t(xx-one*muhat) %*% Ahatinv %*% (xx-one*muhat)
sigmahat <- sqrt(Q0/nn)
#
ML[ss, ] <- c(muhat,sigmahat,rhohat)
print(c(ss,ML[ss, ]))
}}

plot(ML[ ,c(1,3)])

print(round(nn*var(ML),4))
print(round(cor(ML),4))