######################################################################################
#
# This is a modification of the Code developed by Z&W for Chapter 6
#
######################################################################################
#Get started
module(finmetrics)
#Look at the funciton OLS
args(OLS)
######################################################################################
# Example 1: Single index model and CAPM for Microsoft
######################################################################################
# Create the so-called "excess returns" using data
# from the finmetrics files singleIndex.dat and rf.30day
# The first of these holds closing prices of msft and sp500 from Jan 1990 to Jan 2001
# and nominal 30-day T-bill rates.
#First lets look:
colIds(singleIndex.dat)
colIds(rf.30day)
#Now get the returns. In fact, these are only approximate since dividends are ignored.
ret.ts <- getReturns(singleIndex.dat, type="continuous")
#For the short rate we convert the simple returns to the continuous returns
#This does not make much difference, but we want to be compatible with the text.
excessRet.ts <- seriesMerge(ret.ts,log(1+rf.30day))
excessRet.ts[,"MSFT"] <- excessRet.ts[,"MSFT"] - excessRet.ts[,"RF"]
excessRet.ts[,"SP500"] <- excessRet.ts[,"SP500"] - excessRet.ts[,"RF"]
excessRet.ts <- excessRet.ts[,1:2]
#This code makes first makes excessRet.ts a "three column" time series
#then chops it back down to just the two we want.
# Now lets take a look aout our creation
# by plotting the time series returns scatterplot of our "excess returns"
#
par(mfrow=c(2,1))
plot(excessRet.ts, main="Monthly excess returns on Microsoft
and S&P 500 Index",plot.args=list(lty=c(1,3)))
legend(0,-0.2,legend=c("MSFT","S&P 500"),lty=c(1,3))
plot(seriesData(excessRet.ts[,"SP500"]),
seriesData(excessRet.ts[,"MSFT"]),
main="Scatterplot of Returns",
xlab="SP500",ylab="MSFT")
#
# Now we estimate the alpha and beta of of the CAPM ...
#
# First a look at the OLS function
args(OLS)
# Now we specify the modle and create the OLS fit
ols.fit = OLS(MSFT~SP500, data=excessRet.ts)
# Let's look at what we have created...
class(ols.fit)
names(ols.fit)
# Some of these are a little obscure, like R
# Could they mean R-squared?
# No, as it happens. After a little digging, I find that this is the
# R of the "QR method" --- not someting we'll need.
# Still, it reminds us to use the help functions and follow where they lead.
# For example
help(OLS)
# Leads us to OLS.object
help(OLS.object)
# Still no news about R...
# but there is other useful information, like the METHODS
# Let's take a direct look..
ols.fit
# and a look using a "method"
summary(ols.fit)
#
#
# Let's explore the methods...
#
coef(ols.fit)
fitted(ols.fit)[1:3]
resid(ols.fit)[1:3]
vcov(ols.fit)
# We can use these to compute the releavant t-statistics
coef(ols.fit)/sqrt(diag(vcov(ols.fit)))
# This gives us t-values versus a null hypothesis of alpha=0 and beta=0
# Even for mighty Microsoft we don't get a very impressive alpha test
# Still, we can be certain that beta is not zero
# On the other hand ...
# More often we are interested in a null of alpha=0 and beta=1
(coef(ols.fit)-c(0,1))/sqrt(diag(vcov(ols.fit)))
# Pretty solid evidence that beta is not 1
# Due dilligence also insists that we look at the residuals
summaryStats(residuals(ols.fit))
# What's the bottom line here? The mean is "zero" just by the logic of OLS
# The variance is just the variance (easier to understand its square root),
# but there is a message in the skewness and kurtosis. These would be zero
# and 3 in a perfect world. Here skewness is negative and kurtosis is "too large"
# Recall these are for monthly returns. For daily returns, it would be worse.
########################################################
# Advanced section --- to be omitted in 434
#
# Now we could consider the
# construction of Wald tests using the Rb=r framework.
# For 434 this is a bit much so we will skip it,
# but for those who are on top of linear algebra, it is instructive,
# So I'll leave the ZW code...(but comment it out)
#
# Rmat = diag(2)
# rvec = c(0,1)
# bhat = coef(ols.fit)
# avarRbhat = Rmat%*%vcov(ols.fit)%*%t(Rmat)
# wald.stat =
# t(Rmat%*%bhat-rvec)%*%solve(avarRbhat)%*%(Rmat%*%bhat-rvec)
# as.numeric(wald.stat)
# p.value = 1 - pchisq(wald.stat,2)
# p.value
# p.value = 1 - pf(F.stat,2,ols.fit$df.resid)
# p.value
######### End Advanced section #######################
# Suppose we want a JOINT test of the
# hypotheses that alpha=0 and beta=1
# We can use the function waldTest to do this....
waldTest(ols.fit,Intercept==0,SP500==1)
# We can also use the Wald Test to see if a specified
# relationship, or set of relationships, holds for our data
# Here is a one relationship test...
waldTest(ols.fit,Intercept-2*SP500==2)
#
# An OLS object contains information about the likelihood
# This can be used to construct "likelihood ratio tests"
# which can help you decide to use one model over another.
# For the "restricted CAPM" we insist that alpha=0.
# To do this in S we use "-1" in the model as
# language for "leave out the constant term"
ols.fit2 = OLS(MSFT~SP500-1,data=excessRet.ts)
#Now let's look at the log likelihoods for our two fits.
#We use the S function IC with the type="loglike" to extract
#the likelihood from the object ols.fit
IC(ols.fit,"loglike")
IC(ols.fit2,"loglike")
# How can we tell if these are "significantly" different.
# This is where the likelihood ratio test comes in.
# First we define and calculate the LR...
LR = -2*(IC(ols.fit2,"loglike")-IC(ols.fit,"loglike"))
# and look at it...
LR
# Well, duh, this is just a number. What does it mean?
# In this case we know that under the null hypothesis
# LR has a Chi squared distribution with one degree of freedom
# so we can compute a p-value
1 - pchisq(LR,1)
#
# There are several diagnostic plots worth our consideration.
# We get a buffet of these using the plot function for the ols object
plot(ols.fit)
# We can select what we want from the menu, or we
# can select on the code line and prettify it.
plot(ols.fit,which.plot=3)
qqPlot(resid(ols.fit), strip.text="ols.fit",
xlab="Quantile of Standard Normal",
ylab="Residuals",main="Normal QQ Plot")
#
# How about testing the residuals for normality and autocorrelation
#
normalTest(ols.fit,method="jb")
autocorTest(ols.fit,method="lb")
# Looks like we reject normality (as usua) but fail to reject the hypothesis
# "no autocorrelation." We know from experience, that the second "conclusion"
# is lag dependent and the default value is not as sensitive as smaller values.
# Still, it is enough for the moment.
# Could our analyis have somehow depended on the period we chose?
# We can address this by doing a
# subset regression using start= and end= options
#
ols96.fit= OLS(MSFT~SP500,data=excessRet.ts,
start="Jan 1996",in.format="%m %Y")
summary(ols96.fit)
# The coefficients differ, but the story is the same.
# A fancier version of subset regression can be used to just look
# at special days, such as days of positive return on the S&P
olspos.fit=OLS(MSFT~SP500,data=excessRet.ts,
subset=(SP500>=0))
#Zowie...there is some news here. On positive days for S&P MSFT has a smaller beta
#and on negative days it has a larger beta. This sounds like a bad feature.
#Yet for the period in question, MSFT was a fine holding. Mystery?
#
# Estimate January Effect
#
# logical vector picking out January
is.Jan = (months(positions(excessRet.ts))=="Jan")
Jan.dum = timeSeries(pos=positions(excessRet.ts),
data=as.integer(is.Jan))
newdat.ts = seriesMerge(excessRet.ts,Jan.dum)
colIds(newdat.ts)[3] = "Jan.dum"
#
# Plain vanilla dummy variable regression
# to see any excess contribution due to January dummy
summary(OLS(MSFT~SP500+Jan.dum,data=newdat.ts))
# How about adding an interaction term for SP and January?
# we can do this with the colon notation.
OLS(MSFT~SP500+Jan.dum+SP500:Jan.dum,data=newdat.ts)
# or equivalently we can use the "star" notation where A*B means A+B+A:B
summary(OLS(MSFT~SP500*Jan.dum,data=newdat.ts))
# What is the bottom line here? The interaction term looks
# more important that the raw January dummy. How would you interpret this?
# A few lines to illustrate the method "predict"
# for an OLS fit. Here we suppose that a leprechaun tells us
# that the percent return on the SP500 for the next three months
# is -20 percent, zero percent, and +20 percent. How would this
# combine with our model to provide information about MSFT returns?
# First put McFlurry's values into a data frame and give it a name
sp500.new = data.frame(c(-0.2,0,2))
colIds(sp500.new) = "SP500"
# Now use the method "predict" and the OLS object ols.fit to predict the
# next three months (132, 133, and 134) for MSFT
ols.pred = predict(ols.fit,n.predict=3,newdata=sp500.new)
# This object has the goodies, let's check its class...
class(ols.pred)
# and take a default look ...
ols.pred
# and look at an official summary ...
summary(ols.pred)
# and finally a nice picture...
plot(ols.pred,xold=excessRet.ts[,1],n.old=5,width=2)