# Chapter 13 Scripts # Exhibit 13.1 win.graph(width=4.875,height=2.5,pointsize=8) t=1:96 # n=96 cos1=cos(2*pi*t*4/96) cos2=cos(2*pi*(t*14/96+.3)) plot(t,cos1, type='o', ylab='Cosines') lines(t,cos2,lty='dotted',type='o',pch=4) # Exhibit 13.2 # You may want to plot many of the graphs at fullscreen resolution to see more detail win.graph(width=4.875,height=2.5,pointsize=8) t=1:96 # n=96 cos1=cos(2*pi*t*4/96) cos2=cos(2*pi*(t*14/96+.3)) y=2*cos1+3*cos2 plot(t,y,type='o',ylab=expression(y[t])) # Exhibit 13.3 # Assumes that y was created for Exhibit 13.2 sp=periodogram(y) # Note that spec has been altered to default to taper=0 and detrend=F abline(h=0) axis(1,at=c(0.04167,.14583)) # Exhibit 13.4 # Hidden periodicities win.graph(width=4.875,height=2.5,pointsize=8) t=1:96 # n=96 set.seed(134) integer=sample(48,2) # randomly sample 2 integers from 1:48 without replacement freq1=integer[1]/96; freq2=integer[2]/96 A1=rnorm(1,0,2); B1=rnorm(1,0,2) # sample normal "amplitudes" A2=rnorm(1,0,3); B2=rnorm(1,0,3) w=2*pi*t y=A1*cos(w*freq1)+B1*sin(w*freq1)+A2*cos(w*freq2)+B2*sin(w*freq2)+rnorm(96,0,1) plot(t,y,type='o',ylab=expression(y[t])) # Exhibit 13.5 # Assumes that y was created for Exhibit 13.4 above periodogram(y) abline(h=0) # Exhibit 13.6 data(star) plot(star,xlab='Day',ylab='Brightness') # Exhibit 13.7 periodogram(star,ylab='Variable Star Periodogram'); abline(h=0) # Exhibit 13.8 aliasing win.graph(width=4.875, height=2.5,pointsize=8) t=seq(0,8,by=.05) plot(t,cos(2*pi*t/4),axes=F,type='l',ylab=expression(Y[t]), xlab='Discrete Time t') axis(1,at=c(1,2,3,4,5,6,7)) axis(1); axis(2); box() lines(t,cos(2*pi*t*3/4),lty='dashed',type='l') abline(h=0) points(x=c(0:8),y=cos(2*pi*c(0:8)/4),pch=19) # Exhibits 13.9 & 13.10 # MA1 spectral density win.graph(width=4.875, height=2.5,pointsize=8) theta=.9 # Reset theta for other MA(1) plots ARMAspec(model=list(ma=-theta)) # R uses the plus convention in the MA specification. # Exhibits 13.11 # MA(2) models win.graph(width=4.875, height=2.5,pointsize=8) theta1=1; theta2=-0.6 # Reset values of theta1 & theta2 for other MA(2) models ARMAspec(model=list(ma=-c(theta1,theta2))) # Exhibit 13.12 & 13.13 # AR(1) win.graph(width=4.875, height=2.5,pointsize=8) phi=0.9 # Reset value of phi for other AR(1) models ARMAspec(model=list(ar=phi)) # Exhibit 13.14 & 13.15 # AR(2) win.graph(width=4.875, height=2.5,pointsize=8) phi1=1.5; phi2=-.75 # Reset values of phi1 & phi2 for other AR(2) models ARMAspec(model=list(ar=c(phi1,phi2))) # Exhibit 13.17 # ARMA(1,1) win.graph(width=4.875, height=2.5,pointsize=8) theta=0.8 ; phi=0.5 # Reset parameters for other ARMA(1,1) models ARMAspec(model=list(ar=phi,ma=-theta)) # Exhibit 13.18 # Seasonal AR win.graph(width=4.875, height=2.5,pointsize=8) phi=.5 ; PHI=.9 # Reset parameters for other models ARMAspec(model=list(ar=phi,seasonal=list(sar=PHI,peroid=12))) # Exhibit 13.19 # Seasonal MA win.graph(width=4.875, height=2.5,pointsize=8) theta=.4 ; Theta=.9 # Reset parameters for other models ARMAspec(model=list(ma=-theta,seasonal=list(sma=-Theta,period=12))) # Exhibit 13.20 # Simulated AR(1) with phi=-0.6 and n = 200 for Exhibits 13.18,13.26, and 13.27 win.graph(width=4.875,height=2.5,pointsize=8) set.seed(271435) n=200 phi=-0.6 y=arima.sim(model=list(ar=phi),n=n) sp=spec(y,log='no',xlab='Frequency',ylab='Sample Spectral Density',main='', sub='',demean=T) lines(sp\$freq,ARMAspec(model=list(ar=phi),freq=sp\$freq,plot=F)\$spec, lty='dotted') abline(h=0) # Exhibits 13.21,13.22, 13.23, 13.24, 13.25, & 13.26 #Simulation of spectral density sampling properties win.graph(width=4.875,height=2.5,pointsize=8) set.seed(145357) phi1=1.5; phi2=-.75; theta1=0; theta2=0 specm=NULL nrepl=1000 n=96; df=5 # df used only when generating non-normal white noise terms for (i in 1:nrepl){ # Remove the # signs for the next three lines when generating non-normal white noise terms #y=arima.sim(model=list(ar=c(phi1,phi2),ma=c(-theta1,-theta2)), # rand.gen=function(n,...){if(df<=2) stop('df must be greater than 2') else # return(rt(n,...)*sqrt((df-2)/df)) },n=n,df=df) # comment out the line below when generating non-normal white noise terms y=arima.sim(model=list(ar=c(phi1,phi2),ma=c(-theta1,-theta2)),n=n) res=spec(y,plot=F,demean=T) specm=cbind(specm,res\$spec) } dim1=dim(specm)[1] means=apply(specm,1,mean) vars=apply(specm,1,var) cor(t(specm[1:6,])) # Below, compute the empirical frequency of a significant correlation # between two spectral density estimates. mean(abs(cor(t(specm))[outer(1:dim1,1:dim1,"<")])>1.96/sqrt(n)) sdf1=ARMAspec(model=list(ar=c(phi1,phi2),ma=-c(theta1,theta2)),freq=res\$freq, plot=F)\$spec plot(y=sdf1,x=res\$freq, xlab='Frequency',ylab='Average Spectral Density Estimates',type='l', main='',ylim=extendrange(c(means,sdf1))) points(y=means,x=res\$freq) abline(h=0) m=round(n/2) plot(y=sdf1[-m],x=res\$freq[-m],xlab='Frequency', ylab='Standard Deviation of Spectral Density Estimate', type='l',main='',ylim=extendrange(c(sqrt(vars[-m]),sdf1))) points(y=sqrt(vars[-m]),x=res\$freq[-m]) abline(h=0) # Exhibits 13.23 & 13.25 win.graph(width=2.5,height=2.5,pointsize=8) # Replace 40 to obtain any frequency needed qqplot(qchisq(((1:nrepl)-.5)/nrepl,df=2),specm[40,], xlab='Chi-Square Quantitles',ylab='Sample Spectral Distribution')