######################## # Chapter 4 ######################## # The forward equation Different types of behaviour ar2B = function(lambda1,lambda2,n1){ n = (n1+200) gen = rnorm(n) x = rep(0,n) for(j in c(3:n)){ x[j] = ((lambda1+lambda2)*x[j-1] - lambda1*lambda2*x[j-2] + gen[j]) } x1 = x[-c(1:200)] return(x1) } test1 = ar2B(lambda1=0.8,lambda2=0.6,n1=100) test2 = ar2B(lambda1=0.8,lambda2=1.1,n1=100) test3 = ar2B(lambda1=0.8,lambda2=1,n1=100) ######################### # Section 4.3.6 Simulating AR(2) with complex roots. ar2complex = function(lambda,omega,n1){ n = (n1+200) gen = rnorm(n) x = rep(0,n) for(j in c(3:n)){ x[j] = (2*lambda*cos(omega)*x[j-1] - (lambda**2)*x[j-2] + gen[j]) } x1 = x[-c(1:200)] return(x1) } P=pi/3 P=0 test1 = ar2complex(lambda = 0.5,omega = P,n1=200) test2 = ar2complex(lambda = 0.9,omega = P,n1=200) time = c(1:200) par(mfrow=c(1,2)) plot(time,test1,lwd=2,type="l",col="blue",ylim=c(min(test1,test2),max(test1,test2)),ylab="r=0.5") plot(time,test2,lwd=2,type="l",col="red",ylim=c(min(test1,test2),max(test1,test2)),ylab="r=0.9") n=200 F1 = abs(fft(test1)/sqrt(200))**2 F2 = abs(fft(test2)/sqrt(200))**2 F1 = F1[c(1:(n/2))] F2 = F2[c(1:(n/2))] freq = 2*pi*c(1:(n/2))/n par(mfrow=c(1,2)) plot(freq,F1,lwd=2,col="blue",type="l",ylab="Periodogram r = 0.5") lines(c(P,P),c(0,max(F1)),lwd=2,lty=2) plot(freq,F2,lwd=2,col="red",type="l",ylab="Periodogram r = 0.9") lines(c(P,P),c(0,max(F2)),lwd=2,lty=2) ##################################### # Spectral density of AR(2) with complex roots TSpecDenCom2 = function(lambda,omega,n1){ freq = exp(2*pi*1i*c(0:(n1-1))/n1) freq2 = exp(4*pi*1i*c(0:(n1-1))/n1) temp = rep(0,(n1/2)) for(i in c(1:(n1/2))){ temp0 = abs(1-2*lambda*cos(omega)*freq[i]+(lambda**2)*freq2[i]) temp[i] = 1/(temp0**2) } return(temp) } n1=200 P=pi/3 test1 = TSpecDenCom2(lambda = 0.5,omega = P,n1=200) test2 = TSpecDenCom2(lambda = 0.9,omega = P,n1=200) freq = 2*pi*c(0:((n1/2)-1))/n1 plot(freq,test1,type = "l",col="blue",lwd=2,main="spectrum, theta=pi/3, r = 0.5") lines(c(P,P),c(0,max(test1)),lwd=2,lty=2) plot(freq,test2,type = "l",col="red",lwd=2,main="spectrum, theta = pi/3, r = 0.9") lines(c(P,P),c(0,max(test2)),lwd=2,lty=2) P=0 test1 = TSpecDenCom2(lambda = 0.5,omega = P,n1=200) test2 = TSpecDenCom2(lambda = 0.9,omega = P,n1=200) freq = 2*pi*c(0:((n1/2)-1))/n1 plot(freq,test1,type = "l",col="blue",lwd=2,main="spectrum, theta=0, r = 0.5") lines(c(P,P),c(0,max(test1)),lwd=2,lty=2) plot(freq,test2,type = "l",col="red",lwd=2,main="spectrum, theta = 0, r = 0.9") lines(c(P,P),c(0,max(test2)),lwd=2,lty=2) #################################### # Section 4.3.8 ##################################### # SAR models ar12 = function(rho,n1){ n = (n1+1000) gen = rnorm(n) x = rep(0,n) for(j in c(13:n)){ x[j] = (rho*x[j-12] + gen[j]) } x1 = x[-c(1:1000)] return(x1) } # Spectral density for the above model TSpecSAR12 = function(n1,rho){ freq12 = exp(24*pi*1i*c(1:n1)/n1) temp = rep(0,(n1/2)) for(i in c(1:(n1/2))){ temp0 = abs(1-rho*freq12[i]) temp1 = temp0**(-2) # temp2 = abs(1-0.7294*freq[i])**2 temp[i] = temp1 } return(temp) } temp = ar12(rho=0.8,n1=200) time = c(1:200) plot(time, temp,lwd=2,col="blue",type="l") F1 = abs(fft(temp)/sqrt(200))**2 F1 = F1[c(1:(n/2))] spectral = TSpecSAR12(rho=0.8,n1=200) freq = 2*pi*c(1:(n/2))/n par(mfrow=c(1,2)) plot(freq,F1,lwd=1.5,col="blue",type="l",ylab="Periodogram") plot(freq,spectral,lwd=1.5,col="blue",type="l",ylab="spectrum") ############################################## # More advanced SAR model, includes and AR(1) term too ar1.12 = function(phi,rho,n1){ n = (n1+1000) gen = rnorm(n) x = rep(0,n) for(j in c(13:n)){ x[j] = (phi*x[j-1]+rho*x[j-12] + gen[j]) } x1 = x[-c(1:1000)] return(x1) } # Example of another model test = ar1.12(phi=0.3,rho=0.5,m) # Spectral density of the more advanced SAR TSpecSAR1.12 = function(phi,rho,n1){ freq = exp(2*pi*1i*c(1:n1)/n1) freq12 = exp(24*pi*1i*c(1:n1)/n1) temp = rep(0,(n1/2)) for(i in c(1:(n1/2))){ temp0 = abs(1-phi*freq[i]-rho*freq12[i]) temp1 = temp0**(-2) # temp2 = abs(1-0.7294*freq[i])**2 temp[i] = temp1 } return(temp) } ###################################################### # Simulating AR(3), two complex and one real. Using the AR(2) as a basis # Section 4.5 ar3complex = function(lambda,omega,a,n1){ n = (n1+200) # Simulating AR(2) innovations innovations = ar2complex(lambda,omega,n) x = rep(0,n) for(j in c(2:n)){ x[j] = (a*x[j-1] + innovations[j]) } x1 = x[-c(1:200)] return(x1) } # spectral density of above AR(3) TSpecDenAR3 = function(lambda,omega,a,n1){ freq = exp(2*pi*1i*c(0:(n1-1))/n1) freq2 = exp(4*pi*1i*c(0:(n1-1))/n1) temp = rep(0,(n1/2)) for(i in c(1:(n1/2))){ temp0 = abs(1-2*lambda*cos(omega)*freq[i]+(lambda**2)*freq2[i]) temp1 = abs(1-a*freq[i]) temp2 = (temp0**2)*(temp1**2) temp[i] = 1/(temp2) } return(temp) } # Simulating an ARMA(2,2) the AR roots are complex. # b1 and b2 are the MA(2) coefficients arma22complex = function(lambda,omega,b1,b2,n1){ n = (n1+200) gen = rnorm(n) x = rep(0,n) for(j in c(4:n)){ x[j] = (2*lambda*cos(omega)*x[j-1] - (lambda**2)*x[j-2] + gen[j]+b1*gen[j-1]+b2*gen[j-2]) } x1 = x[-c(1:200)] return(x1) } arma32complex = function(lambda,omega,a,b1,b2,n1){ n = (n1+200) # Simulating ARMA(2,2) innovations innovations = arma22complex(lambda,omega,b1,b2,n) x = rep(0,n) for(j in c(2:n)){ x[j] = (a*x[j-1] + innovations[j]) } x1 = x[-c(1:200)] return(x1) } # Spectral density of above model TSpecDenARMA32 = function(lambda,omega,a,b1,b2,n1){ freq = exp(2*pi*1i*c(0:(n1-1))/n1) freq2 = exp(4*pi*1i*c(0:(n1-1))/n1) temp = rep(0,(n1/2)) for(i in c(1:(n1/2))){ temp0 = abs(1-2*lambda*cos(omega)*freq[i]+(lambda**2)*freq2[i]) temp1 = abs(1-a*freq[i]) temp2 = (temp0**2)*(temp1**2) temp3 = abs(1+b1*freq[i]+b2*freq2[i]) temp[i] = (temp3**2)/(temp2) } return(temp) } n1=200 testAR3 = ar3complex(lambda=0.8,omega=pi/3,a=0.6,n1) testARMA32 = arma32complex(lambda=0.8,omega=pi/3,a=0.6,b1=0.5,b2=-0.5,n1) time=c(1:200) PtestAR3 = abs(fft(testAR3)/sqrt(200))**2 PtestARMA32 = abs(fft(testARMA32)/sqrt(200))**2 PtestAR3 = PtestAR3[c(1:(n1/2))] PtestARMA32 = PtestARMA32[c(1:(n1/2))] specAR3 = TSpecDenAR3(lambda=0.8,omega=pi/3,a=0.6,n1) specARMA32 = TSpecDenARMA32(lambda=0.8,omega=pi/3,a=0.6,b1=0.5,b2=-0.5,n1) freq = 2*pi*c(1:(n1/2))/n1 par(mfrow=c(1,2)) plot(time,testAR3,lwd=2,col="blue",type = "l") plot(time,testARMA32,lwd=2,col="red",type = "l") par(mfrow=c(1,2)) plot(freq,PtestAR3,lwd=2,col="blue",type = "l",main = "Periodogram AR(3)") plot(freq,PtestARMA32,lwd=2,col="red",type = "l", main = "Periodogram ARMA(3)") par(mfrow=c(1,2)) plot(freq,specAR3,lwd=2,col="blue",type = "l",main = "Spectral Density AR(3)") plot(freq,specARMA32,lwd=2,col="red",type = "l",main = "Spectral Density ARMA(3,2)") ################################################################# # Simulating an ARFIMA model # load library library("arfima")