?arima.sim
?acf
?adf.test

# Simulating AR(1) processes
# y[t] = rho*y[t-1]+e[t]  
# e[t] ~ N(0,sig^2)
# --------------------------
set.seed(1357)
n   = 10000
sig = 0.1
y   = matrix(0,n,4)
rho = c(0.9,0.99,1.0,1.001)
for (t in 2:n)
    y[t,] = rnorm(4,rho*y[t-1,],sig)

range = range(y[,1:3])
#range = range(y)

par(mfrow=c(1,1))
ts.plot(y[,1],ylim=range,ylab="",lwd=2)
for (i in 2:4)
  lines(y[,i],col=i,lwd=2)
legend("bottomleft",legend=c(paste("rho=",rho[1],sep=""),
paste("rho=",rho[2],sep=""),paste("rho=",rho[3],sep=""),
paste("rho=",rho[4],sep="")),bty="n",col=1:4,lty=1,lwd=2)
abline(h=0,lty=2)

par(mfrow=c(2,2))
for (i in 1:4){
  acf(y[,i],lag.max=50,main="")
  title(paste("rho=",rho[i],sep=""))
}


# Simulating AR(2) processes 
# --------------------------
n   = 100000
sig = 0.1
y   = matrix(0,n,4)

set.seed(654321)
y[,1] = arima.sim(n=n,model=list(ar=c(1.75,-0.85)),sd=sig)
y[,2] = arima.sim(n=n,model=list(ar=c(1.75,-0.7501)),sd=sig)
y[,3] = arima.sim(n=n,model=list(ar=c(0.9,0.0)),sd=sig)
y[,4] = arima.sim(n=n,model=list(ar=c(0.9,0.0999)),sd=sig)
range = range(y)

par(mfrow=c(2,4))
for (i in 1:4) 
  ts.plot(y[,i],ylim=range,ylab="",main=paste("y",i,sep=""))
for (i in 1:4)
  acf(y[,i],lag.max=50,main="")


# Fitting AR(p) for real GDP data
# -------------------------------
#install.packages("tseries")
library("tseries")
data = read.table("https://raw.githubusercontent.com/elara7/R-in-SOE/master/R/R-class/c11/data/q-gdp4708.txt",header=TRUE)

n = nrow(data)
data[c(1,n),]
gdp=log(data[,4])
date = data[,1]+data[,2]/12

par(mfrow=c(2,2))
plot(date,gdp,xlab="Year",ylab="log-gdp",type="l")
plot(date[2:n],diff(gdp),xlab="Year",ylab="diff log-gdp",type="l")
acf(gdp,main="log-gdp")
pacf(diff(gdp),main="diff log-gdp")

m1=ar(diff(gdp),method="mle")

m1 

adf.test(gdp,k=2)


# Fitting AR(p) for real SP500 data
# ---------------------------------
da=read.table("https://raw.githubusercontent.com/elara7/R-in-SOE/master/R/R-class/c11/data/d-sp55008.txt",header=T)
sp5=log(da[,7])
date = da[,1]+(da[,2]*30+da[,3]/31)/366

par(mfrow=c(2,2))
plot(date,da[,7],xlab="Days",ylab="SP500",type="l",main="Price")
plot(date,sp5,xlab="Days",ylab="SP500",type="l",main="Log price")
acf(sp5,main="Log price")
pacf(sp5,main="Log price")

adf.test(sp5,k=10)

m2 = ar(diff(sp5),method="mle")

m2

adf.test(diff(sp5),k=10)