#
#
# 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.
# 
#-----

# fix the random.seed, so I'll always get the same answer:
set.seed(10131985)

nsim = 10^3

# generate random numbers u_k:
U=rnorm(nsim) 

# generate the random numbers x_k:
# 
X = c(0)
for( ii in 1:(nsim-1) ){
  X = c(X,0.8*X[length(X)]+0.6*U[ii])
}


# First plot the sequence U and X:
# 
##postscript("../../WriteUp/Graphics/Chapter2/sect_4_prob_6_raw_data.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,2))
plot(1:nsim,U,col="black")
plot(1:nsim,X,col="black")
##dev.off()

# Next compute the autocorrelations of U and X:
# 
##postscript("../../WriteUp/Graphics/Chapter2/sect_4_prob_6_acf.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,2))
acf(U)
acf(X)
##dev.off()

# The crosscorrelation function:
# 
ccf( U, X )

# Now compute histograms of U and X:
#
##postscript("../../WriteUp/Graphics/Chapter2/sect_4_prob_6_hist.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,2))
hist(U,freq=F,breaks=40)
hist(X,freq=F,breaks=40)
##dev.off()

# Compute the DFT of the u_k and x_k series
#
fft_U = fft(U)
plot( abs(fft_U) )

fft_X = fft(X)
plot( abs(fft_X) )

# Compute the power spectrum of u_k and x_k
# 
aU = acf(U,plot=F)
plot( abs(fft(aU$acf)) )

aX = acf(X,plot=F)
plot( abs(fft(aX$acf)) )

# Compute the cross spectrum of u_k and x_k
# 
ccUX = ccf( U, X, plot=F )
plot( abs(fft(ccUX$acf)) )