# Simulate y as a delayed version of x
# K. Ensor
# STAT 421, Rice University
# ensor@rice.edu
#
# Simulate the process X, from xmodel. Y is then a delayed
#    version of X. The coherence spectrum should be one and
#    the phase spectrum should be a straight line.
#
#

xmodel<-list(ar=c(.5),period=4)
n<-200
temp<-arima.sim(xmodel,n=n)
x<-temp-mean(temp)
ycoef<-c(0,0,.8)  ### For more than 5 lags alter xy below
y<-filter(x,filter=ycoef,method="convolution",sides=1,circular=F)
xy<-cbind(x[6:n],y[6:n])
nxy<-length(xy[,1])

par(mfrow=c(1,1),lwd=3)
tsplot(xy[(nxy-49):nxy,],type="b",col=c(4,3),lty=c(1,6),
	pch=c(3,5)) # Plotting last 50 points, to highlight dependence structure
legend(30,-2,legend=c("X Series","Y Series"),marks=c(3,5),
	col=c(4,3),lty=c(1,6))
acf(xy)
##########################
#Let's take a look at the acf of the prewhitened series

xarfit<-ar(xy[,1],aic=F,order.max=10,method="burg") # AR spectral estimate
yarfit<-ar(xy[,2],aic=F,order.max=10,method="burg")
xypre<-cbind(xarfit$resid[11:nxy],yarfit$resid[11:nxy])
rm(xarfit,yarfit)
acf(xypre)
##########################	

# Plot Estimated Spectral Density

# Frist the AR Spectral Estimate
	par(mfrow=c(2,2),omi=c(.5,1,.5,1),lwd=3)
	bivspecar<-spectrum(xy,method="ar",plot=F)
	matplot(bivspecar$freq,bivspecar$spec, type="l",col=c(4,3),lty=c(1,6))
	title(main="AR Auto Spectral Estimate")

# Second the Nonparametric Smoothed Periodogram Spectral Estimate	  	
	
	bivspecnp<-spectrum(xy,plot=F,method="pgram",spans=15)
	matplot(bivspecnp$freq,bivspecnp$spec, type="l",col=c(4,3),lty=c(1,6))
	title(main="Nonparametric Auto Spectral Estimates")
	legend(-.25,-7,legend=c("X Series","Y Series"),col=c(4,3),lty=c(1,6))

# Now take a look at the estimated squared coherency (both estimates)

	matplot(bivspecnp$freq,cbind((bivspecar$coh)^2,(bivspecnp$coh)^2),
			ylim=c(0,1),type="l",col=c(4,3),lty=c(1,6))
	title(main="Est. Squared Coherencey Spectrum")

# Finally, the estimated phase spectrum (both estimates)

	matplot(bivspecar$freq,cbind(bivspecar$phase,bivspecnp$phase),type="l",
				type="l",col=c(4,3),lty=c(1,6))
	title(main="Est. Phase Spectrum")
	legend(-.35,-.1,legend=c("AR Est","Nonparametric Est"),col=c(4,3),lty=c(1,6))

# What is the slope of the phase spectrum? What is the delay?

	slope<-lsfit(bivspecar$freq,bivspecar$phase, intercept=F)
	slope$coef/(2*pi)
	
	par(mfrow=c(1,1))