setwd("~/teaching/time-series/ts2015spring/data")
###----------------------------------------------------------------------------
### Daily returns
###----------------------------------------------------------------------------

da=read.table("m-gm3dx7508.txt",header=T) # Load data
gm=da[,2]
# Column 2 contains GM stock returns
gm1=ts(gm,frequency=12,start=c(1975,1))
# Creates a ts object.
par(mfcol=c(4,1)) #Put two plots on a page

plot(gm,type='l')
plot(gm1,type='l')
acf(gm,lag=24)
acf(gm1,lag=24)

###----------------------------------------------------------------------------
### CROSS-COVARIANCE AND CORRELATION MATRICES
###----------------------------------------------------------------------------
library("mvtnorm")
library("MTS")
sig=diag(2) # create the 2-by-2 identity matrix
## sig = matrix(c(1, 0.6, 0.6, 1), 2)
x=rmvnorm(300,rep(0,2),sig) # generate random draws
MTSplot(x) # Obtain time series plots (output not shown)
ccm(x)

###----------------------------------------------------------------------------
### Estimation of VAR models.
###----------------------------------------------------------------------------
library("MTS")

da=read.table("q-gdp-ukcaus.txt",header=T)
gdp=log(da[,3:5])
z=gdp[2:126,]-gdp[1:125,]
z=z*100
m1=VAR(z,2)


## Bayesian Estimation
da=read.table("q-gdp-ukcaus.txt",header=T)
x=log(da[,3:5])
dim(x)
dx=x[2:126,]-x[1:125,]
dx=dx*100
C=0.1*diag(7) ### lambda = 0.1
V0=diag(3) ### Vo = I_3
mm=BVAR(dx,p=2,C,V0)


##
da=read.table("m-dec15678-6111.txt",header=T)
head(da)
x=log(da[,2:6]+1)*100  ## compute log returns, in percentages
rtn=cbind(x$dec5,x$dec8)  ## select Decile 5 and 8.
tdx=c(1:612)/12+1961  ## create calendar time
require(MTS) ## loag MTS package
colnames(rtn) <- c("d5","d8")

MTSplot(rtn,tdx)
ccm(rtn,lag=6)

###----------------------------------------------------------------------------
### VMA model
###----------------------------------------------------------------------------
VMAorder(rtn,lag=20)
m1=VMA(rtn,q=1) ## Estimation

MTSdiag(m1)  ### Model checking

names(m1)
r1=m1$residuals
mq(r1,adj=4)

### Exact likelihood estimation
m2=VMAe(rtn,q=1)
MTSdiag(m2)
r2=m2$residuals
mq(r2,adj=4)


### Analysis of monhtly log returns of IBM and KO stocks Sample period: 2001-2011 for 132
### observations Purpose: to demonstrate the difference between exact and conditional
### likelihood estimation in small sample.
da=read.table("m-ibmko-0111.txt",header=T)
head(da)
lrtn=log(da[,2:3]+1)*100
dim(da)
tdx=c(1:132)/12+2001
colnames(lrtn) <- c("ibm","ko")
MTSplot(lrtn,tdx)
ccm(lrtn)
mq(lrtn,10)

yt=diffM(lrtn)
mm=ccm(yt)

m1=VMA(lrtn,q=1,include.mean=F)
m2=VMAe(lrtn,q=1,include.mean=F)
m1=VMA(yt,q=1,include.mean=F)
m2=VMAe(yt,q=1,include.mean=F)

t1=m1$Theta
t2=m2$Theta
eigen(t1)
eigen(t2)

#### Compute the theoretical covariance and CCM matrices
####
phi=matrix(c(.816,-1.116,-.623,1.074,-.643,.615,.592,-.133),2,4)
phi
theta=matrix(c(0,-.801,-1.248,0),2,2)
sig=matrix(c(4,2,2,5),2,2)
VARMAcov(Phi=phi,Theta=theta,Sigma=sig,lag=2)


### Simulation
m1=VARMAsim(400,arlags=c(1,2),malags=c(1),phi=phi,theta=theta,sigma=sig)
names(m1)
zt=m1$series
m2=Eccm(zt,maxp=5,maxq=6)
names(m2)


###----------------------------------------------------------------------------
### VARMA Models
###----------------------------------------------------------------------------
### VARMA modeling
Eccm(zt,maxp=6,maxq=6)
m2=VARMA(zt,p=3,q=1)
m2a=refVARMA(m2,thres=0.8)  ## refine the model
m2b=refVARMA(m2a,thres=1) ## Further refinement
MTSdiag(m2b) ### Model checking