###################################################################################### # # 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)