#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

# Get the data from problem 2.1:
#
# temperature measurements: 
Z = c( 200, 202, 208, 204, 204, 207, 207, 204, 202, 199, 201, 198, 200,
       202, 203, 205, 207, 211, 204, 206, 203, 203, 201, 198, 200, 206,
       207, 206, 200, 203, 203, 200, 200, 195, 202, 204 )

N = length(Z)
q = N/2

a0 = mean(Z)

ts = 1:N 

ai = c()
bi = c()
for( ii in 1:(q-1) ){ 
  fii = ii/N
  ai = c( ai, 2 * sum( cos( 2 * pi * fii * ts ) * Z ) / N )
  bi = c( bi, 2 * sum( sin( 2 * pi * fii * ts ) * Z ) / N )
}

# Fill in the last coefficients:
#
aq = sum( ( (-1)^ts ) * Z )
bq = 0

ai = c( ai, aq )
bi = c( bi, bq )

# Print the analysis of variance table for this timeseries:
#
If = c()
print( sprintf("%4s %10s %8s %15s %8s %12s","ii", "frequency", "period", "periodogram", "D.O.F.", "mean square") )
for( ii in 1:q ){
  fii  = ii/N
  Tii  = 1/fii
  if( ii!=q ){
    Ifii = 0.5 * N * ( ai[ii]^2 + bi[ii]^2 )
    dof = 2
  }else{
    Ifii = N * ai[ii]^2
    dof = 1
  }
  If = c(If, Ifii)
  mean_square = Ifii/dof

  print( sprintf("%4d %10.4f %8.4f %15.4f %8d %12.4f", ii, fii, Tii, Ifii, dof, mean_square) )
}

# Fill in the last part of the table:
#
total             = sum( (Z - mean(Z))^2 )
dof_total         = 2*(q-1) + 1
mean_square_total = total/dof_total 

print( sprintf("                %8s %15.4f %8d %12.4f", "total", total, dof_total, mean_square_total) )

# Lets plot the periodogram:
# 
#postscript("../../WriteUp/Graphics/Chapter2/prob_6_plot_periodogram.eps", onefile=FALSE, horizontal=FALSE)
plot( If )
#dev.off()