Homework3

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

#install.packages("fUnitRoots")
library("fUnitRoots")

## Loading required package: timeDate

## Loading required package: timeSeries

## Loading required package: fBasics

data=read.table("m-unrate-4811.csv", header=T)
gdp=log(data[,1])
head(gdp)

## [1] 1.223775 1.335001 1.386294 1.360977 1.252763 1.280934

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

## [1] 12

adfTest(gdp,lags=12,type=c("c"))

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 12
##   STATISTIC:
##     Dickey-Fuller: -2.8398
##   P VALUE:
##     0.05503 
## 
## Description:
##  Thu Oct 17 20:20:44 2019 by user: Setsnee

dgdp=diff(gdp)
tdx=c(1:length(dgdp))
m2=arima(dgdp,order = c(2,0,0),xreg=tdx)
m2

## 
## Call:
## arima(x = dgdp, order = c(2, 0, 0), xreg = tdx)
## 
## Coefficients:
##          ar1     ar2  intercept  tdx
##       0.1056  0.2374     0.0018    0
## s.e.  0.0352  0.0353     0.0041    0
## 
## sigma^2 estimated as 0.001415:  log likelihood = 1425.87,  aic = -2841.75

m2$coef

##           ar1           ar2     intercept           tdx 
##  1.055739e-01  2.374440e-01  1.803328e-03 -1.307073e-06

sqrt(diag(m2$var.coef))

##          ar1          ar2    intercept          tdx 
## 3.523565e-02 3.526080e-02 4.130150e-03 3.726187e-05

tratio=m2$coef/sqrt(diag(m2$var.coef))
tratio

##         ar1         ar2   intercept         tdx 
##  2.99622398  6.73393746  0.43662524 -0.03507803

vix=log(data$rate)
length(vix)

## [1] 767

#acf(vix)
acf

## function (x, lag.max = NULL, type = c("correlation", "covariance", 
##     "partial"), plot = TRUE, na.action = na.fail, demean = TRUE, 
##     ...) 
## {
##     type <- match.arg(type)
##     if (type == "partial") {
##         m <- match.call()
##         m[[1L]] <- quote(stats::pacf)
##         m$type <- NULL
##         return(eval(m, parent.frame()))
##     }
##     series <- deparse(substitute(x))
##     x <- na.action(as.ts(x))
##     x.freq <- frequency(x)
##     x <- as.matrix(x)
##     if (!is.numeric(x)) 
##         stop("'x' must be numeric")
##     sampleT <- as.integer(nrow(x))
##     nser <- as.integer(ncol(x))
##     if (is.na(sampleT) || is.na(nser)) 
##         stop("'sampleT' and 'nser' must be integer")
##     if (is.null(lag.max)) 
##         lag.max <- floor(10 * (log10(sampleT) - log10(nser)))
##     lag.max <- as.integer(min(lag.max, sampleT - 1L))
##     if (is.na(lag.max) || lag.max < 0) 
##         stop("'lag.max' must be at least 0")
##     if (demean) 
##         x <- sweep(x, 2, colMeans(x, na.rm = TRUE), check.margin = FALSE)
##     lag <- matrix(1, nser, nser)
##     lag[lower.tri(lag)] <- -1
##     acf <- .Call(C_acf, x, lag.max, type == "correlation")
##     lag <- outer(0:lag.max, lag/x.freq)
##     acf.out <- structure(list(acf = acf, type = type, n.used = sampleT, 
##         lag = lag, series = series, snames = colnames(x)), class = "acf")
##     if (plot) {
##         plot.acf(acf.out, ...)
##         invisible(acf.out)
##     }
##     else acf.out
## }
## <bytecode: 0x0000000017199168>
## <environment: namespace:stats>

m1=arima(vix,order=c(0,1,1))
m1

## 
## Call:
## arima(x = vix, order = c(0, 1, 1))
## 
## Coefficients:
##          ma1
##       0.0958
## s.e.  0.0305
## 
## sigma^2 estimated as 0.001509:  log likelihood = 1401.28,  aic = -2798.56

Box.test(m1$residuals,lag=10,type="Ljung")

## 
##  Box-Ljung test
## 
## data:  m1$residuals
## X-squared = 116.25, df = 10, p-value < 2.2e-16

pp=1-pchisq(14.25,9) 
pp

## [1] 0.113706

#plot(vix)

# Exercise 7 
da=read.table("q-jnj-earns-9211.csv",header = T)
head(da)

##   day mon year earns
## 1  30   1 1992  0.11
## 2  23   4 1992  0.18
## 3  21   7 1992  0.18
## 4  20  10 1992  0.17
## 5   1   2 1993  0.12
## 6  29   4 1993  0.20

jnj=da$earns
jnj1=ts(jnj,frequency=12,start=c(1992,1))
#par(mfcol=c(2,1))
#plot(jnj1,xlab="year",ylab="returns")
#title(main="(a): Simple returns")
#acf(jnj,lag=24) # command to obtain sample ACF of the data
ln.jnj=log(jnj+1)
Box.test(jnj,lag=12,type="Ljung")

## 
##  Box-Ljung test
## 
## data:  jnj
## X-squared = 584.59, df = 12, p-value < 2.2e-16

Box.test(ln.jnj,lag=12,type="Ljung")

## 
##  Box-Ljung test
## 
## data:  ln.jnj
## X-squared = 598.48, df = 12, p-value < 2.2e-16

jnj2=ts(jnj,frequency=12,start=c(1992,1))
par(mfcol=c(2,1))
#plot(jnj1,xlab="year",ylab="returns")
#title(main="(a): Simple returns")
#acf(jnj1,lag=24)
gnp=diff(ln.jnj)
dim(da)

## [1] 78  4

tdx=c(1:78)/4+1992
#par(mfcol=c(2,1))
#plot(tdx,ln.jnj,xlab="year",ylab="gnp",type="l")
#plot(tdx[2:78],gnp,type="l",xlab="year",ylab="growth")
#acf(gnp,lag=12)
#pacf(gnp,lag=12) # compute PACF
m1=arima(gnp,order=c(3,0,0))
m1

## 
## Call:
## arima(x = gnp, order = c(3, 0, 0))
## 
## Coefficients:
##           ar1      ar2      ar3  intercept
##       -0.9547  -0.9342  -0.9461     0.0087
## s.e.   0.0324   0.0363   0.0277     0.0006
## 
## sigma^2 estimated as 0.0003623:  log likelihood = 192.15,  aic = -374.3

#tsdiag(m1,gof=12)  # model checking discussed later
p1=c(1,-m1$coef[1:3]) # set-up the polynomial
r1=polyroot(p1) # solve the polynomial equation
r1

## [1]  0.014824+1.019333i -1.017089+0.000000i  0.014824-1.019333i

Mod(r1)

## [1] 1.019441 1.017089 1.019441

k=2*pi/acos(0.014824/1.019441)
k

## [1] 4.037376

mm1=ar(gnp,method="mle")
mm1$order

## [1] 4

names(mm1)

##  [1] "order"        "ar"           "var.pred"     "x.mean"      
##  [5] "aic"          "n.used"       "n.obs"        "order.max"   
##  [9] "partialacf"   "resid"        "method"       "series"      
## [13] "frequency"    "call"         "asy.var.coef"

print(mm1$aic,digits=3)

##       0       1       2       3       4       5       6       7       8 
## 251.685 232.820 226.296  51.073   0.000   1.821   3.812   5.118   0.102 
##       9      10      11      12 
##   1.430   1.319   2.384   4.215

aic=mm1$aic  
length(aic)

## [1] 13

#plot(c(0:12),aic,type="h",xlab="order",ylab="aic")
#lines(0:12,aic,lty=2)

#Forecasting b)
m1=arima(jnj,order=c(0,0,9))
m1

## 
## Call:
## arima(x = jnj, order = c(0, 0, 9))
## 
## Coefficients:
##          ma1     ma2     ma3     ma4     ma5    ma6     ma7     ma8
##       1.0058  1.2457  1.2014  2.1391  1.5133  1.319  1.0087  1.2378
## s.e.  0.1272  0.1732  0.2721  0.2876  0.2363  0.215  0.2427  0.2446
##          ma9  intercept
##       0.3804     0.6094
## s.e.  0.1330     0.0791
## 
## sigma^2 estimated as 0.003708:  log likelihood = 96.24,  aic = -170.47

m1=arima(jnj,order=c(0,0,9),fixed=c(NA,0,NA,0,0,0,0,0,NA,NA))
m1

## 
## Call:
## arima(x = jnj, order = c(0, 0, 9), fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, NA, NA))
## 
## Coefficients:
##          ma1  ma2     ma3  ma4  ma5  ma6  ma7  ma8      ma9  intercept
##       0.1294    0  0.6503    0    0    0    0    0  -0.2064     0.6002
## s.e.  0.2891    0  0.3004    0    0    0    0    0   0.2011     0.0468
## 
## sigma^2 estimated as 0.06913:  log likelihood = -7.43,  aic = 24.87

sqrt(0.06913)

## [1] 0.2629258

Box.test(m1$residuals,lag=12,type="Ljung")

## 
##  Box-Ljung test
## 
## data:  m1$residuals
## X-squared = 394.68, df = 12, p-value < 2.2e-16

pv=1-pchisq(394.68,9)
pv

## [1] 0

pv=1-pchisq(394.68,9)
pv

## [1] 0

m1=arima(jnj[1:78],order=c(0,0,9),fixed=c(NA,0,NA,0,0,0,0,0,NA,NA))
m1

## 
## Call:
## arima(x = jnj[1:78], order = c(0, 0, 9), fixed = c(NA, 0, NA, 0, 0, 0, 0, 0, 
##     NA, NA))
## 
## Coefficients:
##          ma1  ma2     ma3  ma4  ma5  ma6  ma7  ma8      ma9  intercept
##       0.1294    0  0.6503    0    0    0    0    0  -0.2064     0.6002
## s.e.  0.2891    0  0.3004    0    0    0    0    0   0.2011     0.0468
## 
## sigma^2 estimated as 0.06913:  log likelihood = -7.43,  aic = 24.87

predict(m1,10)

## $pred
## Time Series:
## Start = 79 
## End = 88 
## Frequency = 1 
##  [1] 0.9486664 0.6352115 0.8469624 0.5757401 0.4878375 0.5478441 0.4812723
##  [8] 0.5757913 0.4822479 0.6002087
## 
## $se
## Time Series:
## Start = 79 
## End = 88 
## Frequency = 1 
##  [1] 0.2629332 0.2651269 0.2651269 0.3154777 0.3154777 0.3154777 0.3154777
##  [8] 0.3154777 0.3154777 0.3201099

#install.packages("fBasics")
# model comparison
da=read.table("q-gdpc96.csv",header=T)
head(da)

##   year mon day    gnp
## 1 1947   1   1 1780.4
## 2 1947   4   1 1778.1
## 3 1947   7   1 1776.6
## 4 1947  10   1 1804.0
## 5 1948   1   1 1833.4
## 6 1948   4   1 1867.6

#gdp=log(da$gdp)
dgdp=diff(gdp)
m1=ar(dgdp,method='mle')
m1$order

## [1] 12

m2=arima(dgdp,order=c(3,0,0))
m2

## 
## Call:
## arima(x = dgdp, order = c(3, 0, 0))
## 
## Coefficients:
##          ar1     ar2     ar3  intercept
##       0.0674  0.2213  0.1585     0.0013
## s.e.  0.0358  0.0350  0.0359     0.0024
## 
## sigma^2 estimated as 0.001379:  log likelihood = 1435.49,  aic = -2860.99

m3=arima(dgdp,order=c(3,0,0),season=list(order=c(1,0,1),period=4))
m3

## 
## Call:
## arima(x = dgdp, order = c(3, 0, 0), seasonal = list(order = c(1, 0, 1), period = 4))
## 
## Coefficients:
##          ar1     ar2     ar3    sar1     sma1  intercept
##       0.0514  0.2102  0.1604  0.1885  -0.1228     0.0013
## s.e.  0.0370  0.0358  0.0360  0.2306   0.2317     0.0025
## 
## sigma^2 estimated as 0.001374:  log likelihood = 1436.97,  aic = -2859.95