Advanced Statistics with R code 90% doen just check it and make a summary

The average weekly earnings of google shares are $1.005 per share with a standard deviation of $0.045. The distribution of returns is more or less symmetric but have high peak. Normality test (Jarque Bera Test) rejects hypothesis of normality for earnings data at 5% level.

Jarque Bera Test

data: returns

X-squared = 54.1642, df = 2, p-value = 1.731e-12

The distribution of log(returns) is also not closer to normal distribution, with symmetry close to 0 but kurtosis statistics is 1.5 which is far from normal distribution statistic. The Jarque-Bera test is significant (p-value < 0.05) indicating that hypothesis of normal distribution for log(returns) can be rejected.

Jarque Bera Test

data: rets

X-squared = 49.4632, df = 2, p-value = 1.816e-11

The following plot is time series plot of weekly price of google shares which shows increasing trend.

Time series Plot of Google Prices

Time plot of returns show that google stock returns shows periods of high volatility at various times. Most returns are between +/- 10%. Sample moments show that distribution of log returns is somewhat symmetric with very thin tails (kurtosis = 1.56).

To check null hypothesis of non stationarity:

Dickey Fuller Test is expressed as H0: φ1=1 vs Ha: φ1<>1 where the null hypothesis indicates unit-root non-stationarity.

The Dickey-Fuller test shows that the earnings time series is unit-root stationary. The test p-values for lags 7 is less than 0.05, therefore the null hypothesis of non stationarity can be rejected (small p-values< 0.05).

Augmented Dickey-Fuller Test

data: rets

Dickey-Fuller = -7.1264, Lag order = 7, p-value = 0.01

alternative hypothesis: stationary

To check serial correlations in the log returns:

Log returns are not serially correlated as shown by the Ljung-Box test with p-values > 0.05 and the autocorrelation plot.

Box-Pierce test

data: coredata(rets)

X-squared = 0.686, df = 1, p-value = 0.4075

To look for evidence of ARCH effects in the log returns:

The analysis below shows a strong ARCH effect. The squared returns are strongly correlated. The Ljung-Box tests on squared residuals are highly significant with p-values < 0.005, and autocorrelation plots shows large autocorrelations for the first 15 lags.

Box-Pierce test

data: coredata(rets^2)

X-squared = 8.1483, df = 1, p-value = 0.00431

Plot of PACF: There is no correlation till lag 10 which shows there is no AR lag.

To fit an ARMA(0,0)-GARCH(2,1) model for the log returns using a normal- distribution for the error terms:

Fitted Model:

Residual Analysis:

Adjusted Pearson Goodness-of-Fit Test:

------------------------------------

group statistic p-value(g-1)

1 20 28.78 0.06948

2 30 38.57 0.11019

3 40 46.39 0.19380

4 50 51.15 0.38929

Above output shows that error terms are normally distributed as all the P values are greater than 0.05.

To fit ARMA(0,0)-eGARCH(1,1) model with Gaussian distribution:

Fitted Model is:

Residual Analysis:

Adjusted Pearson Goodness-of-Fit Test:

------------------------------------

group statistic p-value(g-1)

1 20 28.95 0.06673

2 30 27.82 0.52739

3 40 46.91 0.18000

4 50 50.72 0.40546

Above output shows that error terms are normally distributed as all the P values are greater than 0.05.

To fit ARMA(0,0)-TGARCH(2,1) model with norm-distribution:

Fitted Model:

Residual Analysis:

Adjusted Pearson Goodness-of-Fit Test:

------------------------------------

group statistic p-value(g-1)

1 20 18.29 0.5030

2 30 31.28 0.3525

3 40 29.17 0.8742

4 50 51.36 0.3813

Above output shows that error terms are normally distributed as all the P values are greater than 0.05.

BEST FIT FOR THE MODEL:

garch21 egarch11 gjrgarch21

Akaike -3.436225 -3.451379 -3.463469

Bayes -3.391975 -3.407129 -3.401519

Shibata -3.436449 -3.451603 -3.463905

Hannan-Quinn -3.418814 -3.433968 -3.439094

From above output it can be said that AR(0)-EGARCH(1,1) and AR(0) – GJR(1,1) both have almost equal (smallest) values for all best fit criterion and hence can be considered as best model.

Forecast Analysis:

RE-FIT MODELS LEAVING 100 OUT-OF-SAMPLE OBSERVATIONS FOR FORECAST and EVALUATE STATISTICS for 5 Step Ahead Forecast.

Model	garch21	egarch11	tgarch21
MSE	0.0009497851	0.0009497741	0.0009499554
MAE	0.0224897000	0.0225260500	0.0225533000
DAC	0.6000000000	0.6000000000	0.6000000000

Since MSE is smallest for GJRGARCH model. We should consider this model for forecasting.

CODE:

library(e1071)

library(tseries)

library(xts)

library(rugarch)

price = scan("clipboard")

pricets=ts(price,frequency=52,start=c(2004,45))

returns=pricets/lag(pricets, -1)

summary(price)

skewness(price)

kurtosis(price)

sd(price)

summary(returns)

skewness(returns)

kurtosis(returns)

sd(returns)

#Normality check

hist(returns, prob=TRUE)

lines(density(returns, adjust=2), type="l")

jarque.bera.test(returns)

#log return time series

rets = log(pricets/lag(pricets, -1))

#Normality check

hist(rets, prob=TRUE)

lines(density(rets, adjust=2), type="l")

jarque.bera.test(rets)

plot.ts(price.ts)

# strip off the dates and just create a simple numeric object (require:xts)

ret = coredata(rets);

# creates time plot of log returns

plot(rets)

summary(rets)

skewness(rets)

kurtosis(rets)

sd(returns)

#ADF test for checking null hypothesis of non stationarity

adf.test(rets)

# Computes Ljung-Box test on returns

Box.test( coredata(rets))

# Computes Ljung-Box test on squared returns to test non-linear independence

Box.test(coredata(rets^2))

# Computes Ljung-Box test on absolute returns to test non-linear independence

Box.test(abs(coredata(rets)))

par(mfrow=c(3,1))

# Plots ACF function of vector data

acf(ret)

# Plot ACF of squared returns to check for ARCH effect

acf(ret^2)

# Plot ACF of absolute returns to check for ARCH effect

acf(abs(ret))

#plot returns, square returns and abs(returns)

par(mfrow=c(3,1))

# Plots ACF function of vector data

plot(rets, type='l')

# Plot ACF of squared returns to check for ARCH effect

plot(rets^2,type='l')

# Plot ACF of absolute returns to check for ARCH effect

plot(abs(rets),type='l')

par(mfrow=c(1,1))

# plots PACF of squared returns to identify order of AR model

pacf(coredata(rets),lag=10)

#specify model using functions in rugarch package

#Fit ARMA(0,0)-GARCH(2,1) model

garch21.spec=ugarchspec(variance.model=list(garchOrder=c(2,1)), mean.model=list(armaOrder=c(0,0)))

#estimate model

garch21.fit=ugarchfit(spec=garch21.spec, data=rets)

garch21.fit

#Fit ARMA(0,0)-eGARCH(1,1) model with Gaussian distribution

egarch11.spec=ugarchspec(variance.model=list(model = "eGARCH",

garchOrder=c(1,1)), mean.model=list(armaOrder=c(0,0)))

#estimate model

egarch11.fit=ugarchfit(spec=egarch11.spec, data=ret)

egarch11.fit

#Fit ARMA(0,0)-TGARCH(2,1) model with norm-distribution

gjrgarch21.spec=ugarchspec(variance.model=list(model = "gjrGARCH",

garchOrder=c(2,1)), mean.model=list(armaOrder=c(0,0)), distribution.model = "norm")

#estimate model

gjrgarch21.fit=ugarchfit(spec=gjrgarch21.spec, data=ret)

gjrgarch21.fit

# compare information criteria

model.list = list(garch21 = garch21.fit,

egarch11 = egarch11.fit,

gjrgarch21 = gjrgarch21.fit)

info.mat = sapply(model.list, infocriteria)

rownames(info.mat) = rownames(infocriteria(garch21.fit))

info.mat

# RE-FIT MODELS LEAVING 100 OUT-OF-SAMPLE OBSERVATIONS FOR FORECAST

# EVALUATION STATISTICS

garch21.fit = ugarchfit(spec=garch21.spec, data=ret, out.sample=100)

egarch11.fit = ugarchfit(egarch11.spec, data=ret, out.sample=100)

tgarch21.fit = ugarchfit(spec=gjrgarch21.spec, data=ret, out.sample=100)

garch21.fcst = ugarchforecast(garch21.fit, n.roll=100, n.ahead=5)

egarch11.fcst = ugarchforecast(egarch11.fit, n.roll=100, n.ahead=5)

tgarch21.fcst = ugarchforecast(tgarch21.fit, n.roll=100, n.ahead=5)

fcst.list = list(garch21=garch21.fcst,

egarch11=egarch11.fcst,

tgarch21=tgarch21.fcst)

fpm.mat = sapply(fcst.list, fpm)

fpm.mat