# tsconcepts.ssc			script file for examples used in Time Series Concepts chapter
#							updated for FM 2.0
#
# created: November 14, 2001
# revised: July 7, 2007
#	added example using new predict method for arima objects
#
####################################################################
# 
# Example: simulate gaussian white noise
#

set.seed(101)				
y = rnorm(100,sd=1)

#
# estimate sample moments
#

y.bar = mean(y)
g.hat = acf(y,lag.max=10,type="covariance",plot=F)
r.hat = acf(y,lag.max=10,type="correlation",plot=F)

#
# plot data and SACF
#

par(mfrow=c(1,2))
tsplot(y,ylab="y")
acf.plot(r.hat)
par(mfrow=c(1,1))

#
# test for normality using normalTest
#

qqnorm(y)
qqline(y)

normalTest(y,method="sw")
normalTest(y,method="jb")

#
# Example: sample moments of daily cc returns on Microsoft stock
#

r.msft = getReturns(DowJones30[,"MSFT"],type="continuous")
r.msft@title = "Daily returns on Microsoft"
sample.2000 = (positions(r.msft) > timeDate("12/31/1999")
& positions(r.msft) < timeDate("1/1/2001"))
par(mfrow=c(2,2))
plot(r.msft[sample.2000],ylab="r.msft")
r.acf = acf(r.msft[sample.2000])
hist(seriesData(r.msft))			
qqnorm(seriesData(r.msft))
par(mfrow=c(1,1))

histPlot(r.msft,strip.text="MSFT monthly return")
qqPlot(r.msft,strip.text="MSFT monthly return")

# test for white noise using autocorTest
args(autocorTest.default)
autocorTest(r.msft,lag.n=10,method="lb")

#
# Example: simulate stationary AR(1) with mu=1, phi=0.75 and sigma=1
# Note: set a seed to create reproducable results
#

set.seed(101)
e = rnorm(100,sd=1)
e.start = rnorm(25,sd=1)
y.ar1 = 1 + arima.sim(model=list(ar=0.75),n=100,
innov=e,start.innov=e.start)
mean(y.ar1)
var(y.ar1)		

#
# compute true ACF
#

gamma.j = rep(0.75,10)^seq(10)

#
# plot data
#

par(mfrow=c(2,2))
tsplot(y.ar1,main="Simulated AR(1)")
abline(h=1)
tsplot(gamma.j,type="h",main="ACF and IRF for AR(1)",
ylab="Autocorrelation",xlab="lag")
tmp = acf(y.ar1,lag.max=10)
par(mfrow=c(1,1))

# compute half life
log(0.5)/log(0.75)

#
# Example: ACF for interest rate spreads (from term structure and exchange rates)
#

uscn.id = 100*(lexrates.dat[,"USCNF"]-lexrates.dat[,"USCNS"])
colIds(uscn.id) = "USCNID"
uscn.id@title = "US/CA 30 day interest rate differential"
par(mfrow=c(2,1))
plot(uscn.id,reference.grid=F)
abline(h=0)
tmp = acf(uscn.id)
par(mfrow=c(1,1))

lag.plot(uscn.id,lags=8,layout=c(2,4))

# 
# example of AR process
#

colIds(varex.ts)
smpl = (positions(varex.ts) > timeDate("12/31/1960"))
irate.real = varex.ts[smpl,"RF.REAL"]
par(mfrow=c(2,2))
acf.plot(acf(irate.real,plot=F))
plot(irate.real,main="Monthly Real Interest Rate")
tmp = acf(irate.real,type="partial")
par(mfrow=c(1,1))

# detrended US industrial production
smpl = (positions(varex.ts) > timeDate("1/1/1960"))
ip = as.matrix(log(seriesData(IP.dat[smpl,])))
trend = 1:nrow(ip)
ipg = residuals(lm(ip~trend))

acf(ipg,type="partial")

# estimate AR model using splus function ar
args(ar)
irate.real = irate.real - mean(irate.real)
ar.fit = ar(irate.real,order.max=6)
names(ar.fit)
ar.fit$order
ar.fit$ar

sum(ar.fit$ar)
abs(polyroot(c(1,-ar.fit$ar)))

# 
# simulate signal + noise model
#

set.seed(112)
eps = rnorm(100,sd=1)
eta = rnorm(100,sd=0.5)
z = cumsum(eta)
y = z + eps
dy = diff(y)
par(mfrow=c(2,2))
tsplot(y,main="Signal plus noise",ylab="y")
tsplot(dy,main="1st difference",ylab="dy")
tmp = acf(dy)
tmp = acf(dy,type="partial")
par(mfrow=c(1,1))

q.val = 0.5^2
rho1 = -1/(q.val + 2)
q.val
rho1
#
# equally weighted portfolio of Dow Jones 30 assets
#

ret.DJ30 = getReturns(DowJones30,type="discrete")
ret.mat = as.matrix(seriesData(ret.DJ30))
wts = rep(1/30,30)
ret.port = ret.mat%*%wts

#
# create overlapping annual returns for msft
#

sp500.mret = getReturns(singleIndex.dat[,"SP500"],
type="continuous")
sp500.mret@title = "Monthly returns on S&P 500 Index"
sp500.aret = aggregateSeries(sp500.mret,moving=12,FUN=sum)
sp500.aret@title = "Monthly Annual returns on S&P 500 Index"
par(mfrow=c(2,2))
plot(sp500.mret)
plot(sp500.aret)
acf.plot(acf(sp500.aret,plot=F))
acf.plot(acf(sp500.aret,type="partial",plot=F))
par(mfrow=c(1,1))

#
# estimate ARMA(1,1) model for uscn.id
#

uscn.id.dm = uscn.id - mean(uscn.id)
arma11.mod = list(ar=0.75,ma=0)
arma11.fit = arima.mle(uscn.id.dm,model=arma11.mod)
class(arma11.fit)
names(arma11.fit)
arma11.fit
std.errs = sqrt(diag(arma11.fit$var.coef))
names(std.errs) = colIds(arma11.fit$var.coef)
std.errs

# EZ: 07/05/2007
# note: in splus 8, there is a summary method for objects of class
# "arima" which shows standard errors and t-stats. There is also
# a residuals() function
summary(arma11.fit)
plot(residuals(arma11.fit), reference.grid=F, 
     main="Residuals from ARMA(1,1) fit")
abline(h=0)
seriesPlot(residuals(arma11.fit), one.plot=T, ref.line=T)

#
# estimate ARMA(1,1) with a constant
#

arma11.fit2 = arima.mle(uscn.id,model=arma11.mod,xreg=1)
arma11.fit2
std.errs = sqrt(diag(arma11.fit$var.coef))
names(std.errs) = colIds(arma11.fit$var.coef)
std.errs

plot(arma11.fit)
graphsheet()

# produce h-step ahead forecasts
end(uscn.id.dm)
fcst.dates = timeSeq("7/1/1996","6/1/1997",
by="months",format="%b %Y")
uscn.id.dm.fcst = arima.forecast(uscn.id.dm,n=12,
model=arma11.fit$model,future.positions=fcst.dates)
names(uscn.id.dm.fcst)
uscn.id.dm.fcst[[1]]

# EZ: 7/5/2007. Problem in arima.forecast when timeSeries is 
# extended to accomodate forecast observations. Need to coerce
# NA data to class matrix in order for function to work

	n <- length(fcst.dates)
		x.future <- timeSeries(data = as.matrix(rep(NA, n)), 
			positions = fcst.dates)
		x.fore <- concat(uscn.id.dm, x.future)

# EZ: 7/5/2007
# In FM 3.0, there is a predict() method for objects of class
# "arima" which produces a FM object of class "forecast"
args(predict.arima)
uscn.id.dm.fcst = predict(arma11.fit, n.predict=12)
class(uscn.id.dm.fcst)
uscn.id.dm.fcst
summary(uscn.id.dm.fcst)
plot(uscn.id.dm.fcst, xold=uscn.id.dm, n.old=6)



smpl = positions(uscn.id.dm) >= timeDate("6/1/1995")
plot(uscn.id.dm[smpl,],uscn.id.dm.fcst$mean,
uscn.id.dm.fcst$mean+2*uscn.id.dm.fcst$std.err,
uscn.id.dm.fcst$mean-2*uscn.id.dm.fcst$std.err,
plot.args=list(lty=c(1,2,3,3)))

# 
# estimate AR(2) using OLS
#

ar2.fit = OLS(USCNID~ar(2),data=uscn.id)
ar2.fit
abs(polyroot(c(1,-ar2.fit$coef[2:3])))

# simulate ARCH(1) process

rt = simulate.garch(model=list(a.value=0.1, arch=0.8), 
     n=250, rseed=196)
class(rt)
names(rt)
par(mfrow=c(2,1))
tsplot(rt$et,main="Simulated returns")
tsplot(rt$sigma.t,main="Simulated volatility")
par(mfrow=c(1,1))

#
# estimating the long-run variance of a time series
# y(t) = 1 + 0.75y(t-1) + e(t)

lrv.true = 1/((1-0.75)^2)
set.seed(101)
e = rnorm(100,sd=1)
y.ar1 = 1 + arima.sim(model=list(ar=0.75),innov=e)

# autoregressive estimate
ar1.fit = OLS(y.ar1~ar(1))
rho.hat = coef(ar1.fit)[2]
sig2.hat = sum(residuals(ar1.fit)^2)/ar1.fit$df.resid
lrv.ar1 = sig2.hat/(1-rho.hat)^2
as.numeric(lrv.ar1)
se.ybar.ar1 = sqrt(lrv.ar1/100)

# Newey-West estimate
args(asymp.var)
4*(100/100)^(2/9)
lrv.nw = asymp.var(y.ar1,bandwidth=4)
lrv.nw
sqrt(lrv.nw/100)

#
# compute long-run variance for overlapping returns
#

# simulate random walk
set.seed(201)
e = rnorm(100)
z.rw = cumsum(e)

#
# Example: Variance ratios
#
args(varRatioTest)
# note: x is log prices

# Variance ratios for Dow Jones Industrial average from 1960 - 1990
# note: if hetero=F then first variance ratio is significant
VR.djia = varRatioTest(djia[timeEvent("1/1/1960","1/1/1990"),"close"],
n.periods=60,unbiased=T,hetero=T)
class(VR.djia)
names(VR.djia)
VR.djia
plot(VR.djia)

# Variance ratios for individual Dow Jones stocks from 1991 - 2001
VR.DJ30 = varRatioTest(DowJones30,n.periods=5,unbiased=T,hetero=T)
plot(VR.DJ30)



#
# Example: simulate trend stationary AR(1) with mu=1,delta=0.25, phi=0.75 and sigma=1
#

set.seed(101)
y.tsar1 = 1 + 0.25*seq(100) + 
arima.sim(model=list(ar=0.75),n=100)
tsplot(y.tsar1,ylab="y")
abline(a=1,b=0.25)

set.seed(101)
u.ar1 = arima.sim(model=list(ar=0.75),n=100)
y1 = cumsum(u.ar1)
y1.d = 1 + 0.25*seq(100)+ y1
y2 = rep(0,100)
for (i in 3:100) {
	y2[i] = 2*y2[i-1] - y2[i-2] + u.ar1[i]
}
par(mfrow=c(2,2))
tsplot(u.ar1,main="I(0) innovations")
tsplot(y1,main="I(1) process")
tsplot(y1.d,main="I(1) process with drift")
abline(a=1,b=0.25)
tsplot(y2,main="I(2) process")
par(mfrow=c(1,1))

#
# plot financial data
#

uscn.spot = lexrates.dat[,"USCNS"]
sp500 = singleIndex.dat[,"SP500"]
par(mfrow=c(2,2))
plot.timeSeries(uscn.spot,main="Log US/CA spot exchange rate",reference.grid=F)
plot.timeSeries(log(sp500),main="Log S&P 500 index",reference.grid=F)
plot.timeSeries(rf.30day,main="Nominal 30 day T-bill rate",reference.grid=F)
plot.timeSeries(log(CPI.dat),main="Log of US CPI",reference.grid=F)
par(mfrow=c(1,1))

#
# long memory example
#

# simulate fractional white noise process
set.seed(394)
y.fwn = simulate.FARIMA(list(d=0.3), 500)
par(mfrow=c(2,1))
tsplot(y.fwn)
tmp = acf(y.fwn,lag.max=50)
par(mfrow=c(1,1))

msft.aret = abs(getReturns(DowJones30[,"MSFT"]))
uscn.id = 100*(lexrates.dat[,"USCNF"]-
lexrates.dat[,"USCNS"])
colIds(uscn.id) = "USCNID"

par(mfrow=c(2,1))
tmp = acf(msft.aret,lag.max=100)
tmp = acf(uscn.id,lag.max=50)
par(mfrow=c(1,1))

#########################################################################################
# 								multivariate time series
#########################################################################################

# 
# construct time series of first four asset returns in DowJones30.dat
#
Y = getReturns(DowJones30[,1:4],type="continuous")
colIds(Y)

# compute means and variances and correlations
colMeans(seriesData(Y))			# colMeans doesn't have method for timeSeries
var(Y)								# var does have a method for timeSeries!
cor(Y)
colVars(seriesData(Y))
colStdevs(seriesData(Y))

#
# compute autocovariance and autocorrelation matrices
#

Ghat = acf(Y[,1:2],lag.max=5,type="covariance",plot=F)
Rhat = acf(Y[,1:2],lag.max=5,plot=F)
Rhat
acf.plot(Rhat)
Ghat$acf[1,,]
Rhat$acf[1,,]
Ghat$acf[2,,]
Rhat$acf[2,,]
par(mfrow=c(1,1))

# 
# simulate VAR(1) model: note can't use filter() for efficiency
#

nobs = 250
e.var = rmvnorm(nobs,sd=c(1,1),rho=0.5)
e.var = t(e.var)
pi1 = matrix(c(0.7,0.2,0.2,0.7),2,2,byrow=T)		# VAR(1) coef matrix
mu.vec = c(1,2)											# mean vector 
c.vec = as.vector((diag(2)-pi1)%*%mu.vec) 			# intercept in VAR
y.var = matrix(0,2,nobs)
y.var[,1] = mu.vec										# initial value is mean
for (i in 2:nobs) {
	y.var[,i] = c.vec + pi1%*%y.var[,i-1]+e.var[,i]
}
y.var = t(y.var)
dimnames(y.var) = list(NULL,c("y1","y2"))
eigen(pi1,only.values=T)
colMeans(y.var)
tsplot(y.var)

#
# compute long-run variance
#

M.T = floor(4*(nrow(Y)/100)^(2/9))
lrv.nw = asymp.var(Y,bandwidth=M.T)
lrv.nw
sqrt(diag(lrv.nw))