# Simulation of AR(p) processes. Condition on Phi for causality
# new functions: polyroot to find the complex roots of a polynomial
#                ar.yw to fit an AR model to data (will be explained later)

help(polyroot)
help(ar.yw)

# random choice of AR(2) coefficients in the causality region

u <- runif(2)
phi <- u[1]*c(0,1) + (1-u[1])*(u[2]*c(-2,-1) + (1-u[2])*c(2,-1))
round(phi,4)
z <- polyroot(c(1,-phi))
z
Mod(z)
Arg(z)/(2*pi)
plot(arima.sim(n=400,model=list(ar=phi)))

# AR(2) coefficients from the roots of the polynomial Phi
r <- 1.05
nu <- 1/11
phi <- c(2*cos(2*pi*nu)/r,-1/r^2)
phi
z <- polyroot(c(1,-phi))
z
Mod(z)
Arg(z)/(2*pi)
plot(arima.sim(n=400,model=list(ar=phi)))

# AR coefficients chosen to fit real data
lynx.ar <- ar.yw(log(lynx))
lynx.coef <- lynx.ar$ar
lynx.coef
sort(Arg(polyroot(c(1,-lynx.coef)))/(2*pi))
sort(Mod(polyroot(c(1,-lynx.coef))))
plot(polyroot(c(1,-lynx.coef)))
lines(cos(0.02*pi*(0:100)),sin(0.02*pi*(0:100)))
par(mfrow=c(2,1))
plot(log(lynx))
plot(arima.sim(n=length(lynx),model=list(ar=phi)))
par(mfrow=c(1,1)
    
# Verify that U_t = X_t - \sum_{j=1}^p  phi_j X_t-j (1)
# and X_t = \sum_{j=1}^\infty \psi_j U_{t-j} (2)
# are inverse relations
p <- length(phi)
polyroot(c(1,-phi))
psi <- c(rep(0,p-1),rep(1,501))  #computing psi
for (i in ((p+1):(500+p))) psi[i] <- sum(phi*psi[(i-1):(i-p)])
plot(0:500,psi[p:(500+p)],type="h")
u <- rnorm(700)
#AR(p)
k <- 100  #truncation point for the infinite series in (2)
z <- filter(u,psi[p:(p+k)],"convolution",sides=1) # causal representation (2)
plot(z)
v <- filter(z,c(1,-phi),"convolution",sides=1) # computing innovations
plot(v-u) #should be zero
max(abs(v-u)[(700-k-2):700])
#MA(p)
k <- 100
z <- filter(u,c(1,-phi),"convolution",sides=1) # definition of MA(p) process
plot(z)
v <- filter(z,psi[p:(p+k)],"convolution",sides=1) # Inversion (2)
plot(v-u) #should be zero
max(abs(v-u)[(700-k-2):700])