########################

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