#  Stochastic Calculus, Courant Institute, NYU, Fall 2014
#  http://www.math.nyu.edu/faculty/goodman/teaching/StochCalc2014/index.html

#  Assignment 2, Gaussian processes, discrete time, linear

#  (your name and email here)

printM <- function( M){   # print an m X n matrix}
   dimensions = dim(M)
   m = dimensions[1]      # number of rows
   n = dimensions[2]      # number of columns
   
   for ( i in (1:m)){
      cat( sprintf( 'row %3d: ', i))
      for ( j in (1:n)){
         cat( sprintf(' %6.2f ',M[i,j]) )
        }
    cat('\n')
  }
 }

cat("Welcome to R world\n")

d = 10
T = 200
A = matrix(0., d, d)  # create a d X d matrix of all zeros

alpha = -.35   # above the diagonal
beta  =  .5    # the diagonals
gamma = -.15   # below the diagonal

A[1,1] = beta
A[1,2] = alpha
for ( i in (2:d-1)){
   A[i,i-1] = gamma
   A[i,i]   = beta
   A[i,i+1] = alpha
  }
A[d,d-1] = gamma
A[d,d]   = beta

R = matrix(0., d, d)

r = eigen( A, symmetric = FALSE, only.values = TRUE)  # returns a data structure
lam = r$values                                        # dig out the eigenvalues
spectral.gap = 1. - max( Mod(lam))
cat( sprintf( 'spectral gap = %8.2e\n', spectral.gap))
cat('Eigenvalues: ')
for ( i in (1:d)){
   cat( sprintf( ' %6.3f ', lam[i]))
  }
cat('\n')


C = matrix(0., d, d)
for ( i in (1:d)){
   C[i,i] = 1.
   R[i,i] = 1.
  }

cat('Matrix A is \n')
printM(A)

cat('Matrix R is \n')
printM(R)

cat('Starting covariance is \n')
printM(C)

Cn = array( NA, c(d, d, T))

rms.old = 1.

Cn[,,1] = C
for ( n in (2:T)){
   C.old   = Cn[,,n-1]
   C.new   = A %*% C.old %*% t(A) + R
   Cn[,,n] = C.new
   C.diff  = C.new - C.old
   rms.new = sqrt( sum( diag( ( C.diff%*%t(C.diff)))))
   convergence.rate = (rms.old - rms.new)/rms.old
   cat( sprintf( 'step %4d, ',                             n))
   cat( sprintf( 'rms of covariance differences = %8.3e, ',rms.new))
   cat( sprintf( 'convergence rate = %8.3e\n',             convergence.rate))
   rms.old = rms.new
#   cat( sprintf( 'Covariance after %3d steps is \n', n))
#   printM(C.new)
  }

mp = 2.*spectral.gap - spectral.gap^2
cat( sprintf( 'Spectral gap predicted convergence rate = %8.3e\n', mp))

L = chol(R)

N = 1000  # number of paths to simulate
C.hat = matrix( 0, d, d)

for ( path in (1:N)){

   X = rnorm(d)            # standard normal starting vector
   for ( n in (2:T)){
      Z = rnorm(d)
      X = A%*%X + L%*%Z
     }
   
   C.hat = C.hat + X%*%t(X)
  }

C.hat = C.hat/N

cat( 'Theoretical covariance matrix\n')
printM(C.new)
cat( 'Empirical covariance matrix\n')
printM(C.hat)