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