## An example of Indented code:
M = array(0, c(5,5,5) )
for( i in 1:5 )
{
  for( j in 1:5 )
  {
    for( k in 1:5) 
    {
      M[i,j,k]=i+j+k
    }
  }
}

## An example of code NOT indented
M = array(0, c(5,5,5) )
for( i in 1:5 )
{
for( j in 1:5 )
{
for( k in 1:5) 
{
M[i,j,k]=i+j+k
}
}
}









reps = 1e5
theta = 0.5
p=8

## Generate reps realizations of Z_0, ..., Z_p
Z = array(rnorm(reps*(p+1)),c(reps,p+1))

## Find X in terms of Z
X = array(0, c(reps,p))
for (i in 1:p)
{
  ## X_i = theta * Z_i + Z_{i-1}  (where the Z is one column off)
  X[,i] = theta*Z[,i+1] + Z[,i]
}

## Compute the sample covariance matrix 
SigmaHat = cov(X)
round(SigmaHat,3)






reps = 1e5
theta = 0.5
p=8

## Generate reps realizations of Z_0, ..., Z_p
Z = array(rnorm(reps*(p+1)),c(reps,p+1))

## Find X in terms of Z
X = array(0, c(reps,p))
for (i in 1:p)
{
  ## X_i = theta * Z_i + Z_{i-1}  (where the Z is one column off)
  X[,i] = theta*Z[,i+1] + Z[,i]
}

## Compute the sample covariance matrix 
SigmaHat = cov(X)

## Find the simulated var(Xbar)
Xbar = apply(X, 1, mean)
var(Xbar)

## Find the exact var(Xbar)
exact = (1/p)*(theta^2 + 2*theta +1) - 2*theta*p^(-2)
exact







n=20
reps = 1e5
theta = 5
X = array(runif(n*reps)*theta,c(n,reps))
theta.hat = (n+1)/n * apply(X, 2, max)
sim.var = var(theta.hat)
sim.var

formula.var = theta^2 * ((n+1)^2 / (n*(n+2)) - 1)
formula.var





nrep = 1e4
n = 20
theta = c(2,5,10)

CP = NULL
for (k in 1:3) {
    X = runif(nrep*n, min=0, max=theta[k])
    X = array(X, c(nrep,n))
    M = apply(X, 1, max)
    theta.hat = (n+1)*M/n
    LB = theta.hat - 1.96*sqrt(theta.hat^2 *( (n+1)^2/(n*(n+2)) -1 ) )
    UB = theta.hat + 1.96*sqrt(theta.hat^2 *( (n+1)^2/(n*(n+2)) -1 ) )

    CP[k] = mean( (LB<theta[k]) & (theta[k]<UB) )
}