[Eventstudies-commits] r162 - / pkg pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 6 13:50:40 CET 2014
Author: vikram
Date: 2014-02-06 13:50:39 +0100 (Thu, 06 Feb 2014)
New Revision: 162
Added:
pkg/R/lmAmm.R
Modified:
pkg/NAMESPACE
pkg/R/AMM.R
todo.org
Log:
Added new lmAMM function; fixing the AMM function then moving on to event studies; Work in progress
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2014-01-22 10:01:36 UTC (rev 161)
+++ pkg/NAMESPACE 2014-02-06 12:50:39 UTC (rev 162)
@@ -2,8 +2,14 @@
remap.cumsum, remap.cumprod, remap.event.reindex, ees, eesPlot)
export(marketResidual,
excessReturn)
-export(AMM,
- onefirmAMM,
- manyfirmsAMM,
- firmExposures,
- makeX)
+export(manyfirmsAMM, lmAMM, makeX)
+
+
+S3method(summary, amm)
+S3method(summary, es)
+
+S3method(print, amm)
+S3method(print, es)
+
+S3method(plot, amm)
+S3method(plot, es)
Modified: pkg/R/AMM.R
===================================================================
--- pkg/R/AMM.R 2014-01-22 10:01:36 UTC (rev 161)
+++ pkg/R/AMM.R 2014-02-06 12:50:39 UTC (rev 162)
@@ -184,7 +184,7 @@
exposures <- matrix(NA,nrow=nfirms,ncol=nperiods*ncol(regressors))
rownames(exposures) <- colnames(regressand)
- tmp <- NULL
+ tmp <- NULL
for(i in 1:length(periodnames)){
for(j in 1:ncol(regressors)){
tmp <- c(tmp, paste(colnames(regressors)[j],
Added: pkg/R/lmAmm.R
===================================================================
--- pkg/R/lmAmm.R (rev 0)
+++ pkg/R/lmAmm.R 2014-02-06 12:50:39 UTC (rev 162)
@@ -0,0 +1,300 @@
+
+
+########################
+# Many firms AMM
+########################
+manyfirmsAMM <-
+function(regressand,regressors,
+ lags,dates=NULL, periodnames=NULL,verbose=FALSE){
+ require("doMC")
+ registerDoMC()
+ if(is.null(dates)){
+ dates=c(start(regressors),end(regressors))
+ periodnames="Full"
+ }
+ nperiods <- length(periodnames)
+ if(length(dates) != (nperiods+1)){
+ cat("Mistake in length of dates versus length of periods.\n")
+ return(NULL)
+ }
+ nfirms <- ncol(regressand)
+
+ # Let's get "exposure' and 'sds'. Setting up structures:-
+
+ exposures <- matrix(NA,nrow=nfirms,ncol=nperiods*ncol(regressors))
+ rownames(exposures) <- colnames(regressand)
+ tmp <- NULL
+ for(i in 1:length(periodnames)){
+ for(j in 1:ncol(regressors)){
+ tmp <- c(tmp, paste(colnames(regressors)[j],
+ periodnames[i],sep="."))
+ }
+ }
+ colnames(exposures) <- tmp
+ sds <- exposures
+ colnames(sds) <- paste("sd",colnames(exposures),sep=".")
+
+ # Setup a list structure for an OLS that failed
+ empty <- list(exposures=rep(NA,ncol(regressors)),
+ s.exposures=rep(NA,ncol(regressors)))
+
+ for(i in 1:ncol(regressand)){
+ cat("Doing",colnames(regressand)[i],"\n")
+ if (verbose) {cat ("Doing", colnames(regressand)[i], "\n")}
+ firm.returns <- regressand[,i]
+ dataset <- cbind(firm.returns, regressors) # This is the full time-series
+ this.exp <- this.sds <- NULL
+ for(j in 1:nperiods){ # now we chop it up
+ t1 <- dates[j]
+ t2 <- dates[j+1]
+ this <- window(dataset,start=t1, end=t2)
+ fe <- lmAMM(this[,1],this[,-1],nlags=lags,verbose)
+ if(is.null(fe)) {fe <- empty}
+ this.exp <- c(this.exp, fe$exposures)
+ this.sds <- c(this.sds, fe$s.exposures)
+ }
+ exposures[colnames(regressand)[i],] <- this.exp
+ sds[colnames(regressand)[i],] <- this.sds
+ }
+ list(exposures=exposures, sds=sds, sig=exposures/sds)
+}
+
+###############################################
+## Estimating one firm's exposure in one period
+###############################################
+lmAMM <- function(firm.returns, X, nlags=NA, verbose=FALSE) {
+ do.ols <- function(nlags) {
+ tmp <- cbind(firm.returns, X[,1]) # Assume 1st is stock index, and no lags are required there.
+ labels <- c("firm.returns","market.returns")
+ if (NCOL(X) > 1) {
+ for (i in 2:NCOL(X)) {
+ for (j in 0:nlags) {
+ tmp <- cbind(tmp, lag(X[,i], -j))
+ labels <- c(labels, paste(colnames(X)[i], j, sep="."))
+ }
+ }
+ }
+ tmp <- na.omit(tmp)
+ if (nrow(tmp) < 30) { # refuse to do the work.
+ return(NULL) # returns out of do.ols() only
+ }
+
+ colnames(tmp) <- labels # So the OLS results will look nice
+ lm(firm.returns ~ ., data=as.data.frame(tmp))
+ }
+
+ if (is.na(nlags)) {
+ if (verbose) {cat("Trying to find the best lag structure...\n")}
+ bestlag <- 0
+ bestm <- NULL
+ bestAIC <- Inf
+ for (trylag in 0:min(10,log10(length(firm.returns)))) {
+ thism <- do.ols(trylag)
+ thisAIC <- AIC(thism, k=log(length(thism$fitted.values)))
+ if (verbose) {cat(trylag, " lags, SBC = ", thisAIC, "\n")}
+ if (thisAIC < bestAIC) {
+ bestlag <- trylag
+ bestAIC <- thisAIC
+ bestm <- thism
+ }
+ }
+ nlags <- bestlag
+ m <- bestm
+ } else {
+ m <- do.ols(nlags)
+ if (is.null(m)) {return(NULL)}
+ }
+ # In either event, you endup holding an "m" here.
+ if (verbose) {cat("\n\nThe OLS:\n"); print(summary(m))}
+
+ # Compute a series of exposure measures, and their standard errors.
+ beta <- m$coefficients
+ Sigma <- vcovHAC(m)
+ # First the market.returns
+ exposures <- beta[2] # no lags for market.returns
+ s.exposures <- sqrt(Sigma[2,2])
+ # From here on, there's a block of 1+nlags coeffs for each
+ # of the non-market.returns regressors.
+ if (NCOL(X) > 1) {
+ for (i in 2:NCOL(X)) {
+ n.block1 <- 2 + ((i-2)*(1+nlags)) # Just 2 for the 1st case.
+ n.block2 <- length(beta) - n.block1 - (1 + nlags)
+ w <- c(rep(0, n.block1), rep(1, 1+nlags), rep(0, n.block2))
+ exposures <- c(exposures, w %*% beta)
+ s.exposures <- c(s.exposures, sqrt(w %*% Sigma %*% w))
+ }
+ }
+ results <- m
+ names(exposures) <- names(s.exposures) <- colnames(X)
+ results$exposures <- exposures
+ results$s.exposures <- s.exposures
+ results$nlags <- nlags
+ class(results) <- "amm"
+ return(results)
+}
+
+
+###########################
+# Maintaing NAs in AR model
+###########################
+ARinnovations <- function(x) {
+ stopifnot(NCOL(x) == 1)
+ dt <- NULL
+ if (class(x) == "zoo") {
+ dt <- index(x)
+ x <- as.numeric(x)
+ }
+ non.na.locations <- !is.na(x)
+ x <- x[non.na.locations]
+ m <- ar(x)
+ e <- m$resid
+ # Relocate with the old NA structure
+ result <- rep(NA, length(non.na.locations))
+ result[non.na.locations] <- e
+ if (!is.null(dt)) {result <- zoo(result, order.by=dt)}
+ list(result=result, m=m)
+}
+
+# ---------------------------------------------------------------------------
+# The workhorse called by makeX to return a nice matrix of RHS
+# variables to be used in an analysis.
+do.one.piece <- function(market.returns, others, switch.to.innov, market.returns.purge, nlags, verbose=FALSE) {
+ thedates <- index(market.returns)
+ if (verbose) {
+ cat(" Doing work for period from ",
+ as.character(head(thedates,1)), " to ",
+ as.character(tail(thedates,1)), "\n")
+ }
+ otherlags <- NULL
+ for (i in 1:NCOL(others)) {
+ if (switch.to.innov[i]) {
+ a <- ARinnovations(others[,i])
+ innovated <- a$result
+ if (verbose) {
+ cat(" AR model for ", colnames(others, do.NULL=FALSE)[i], "\n")
+ print(a$m)
+ }
+ otherlags <- c(otherlags, a$m$order)
+ } else {
+ innovated <- others[,i]
+ otherlags <- c(otherlags, 0)
+ }
+ if (i > 1) {
+ innov <- cbind(innov, innovated)
+ } else {
+ innov <- innovated
+ }
+ }
+ if (NCOL(innov) > 1) {colnames(innov) <- colnames(others)}
+ market.returns.purged <- market.returns
+ if (market.returns.purge) {
+ firstpass <- TRUE
+ for (i in 1:NCOL(innov)) {
+ for (j in 0:nlags) {
+ if (firstpass) {
+ z <- lag(innov[,i],-j)
+ labels <- paste(colnames(innov,do.NULL=FALSE)[i],j,sep=".")
+ firstpass <- FALSE
+ } else {
+ z <- cbind(z, lag(innov[,i], -j))
+ labels <- c(labels, paste(colnames(innov,do.NULL=FALSE)[i],j,sep="."))
+ }
+ }
+ }
+ if (NCOL(z) > 1) {colnames(z) <- labels}
+ m <- lm(market.returns ~ ., as.data.frame(cbind(market.returns, z)))
+ if (verbose) {
+ cat(" Model explaining market.returns:\n")
+ print(summary(m))
+ }
+ how.many.NAs <- nlags + max(otherlags)
+ market.returns.purged <- zoo(c(rep(NA,how.many.NAs),m$residuals),
+ order.by=thedates)
+ }
+ # if (verbose) {cat(" Finished do.one.piece()\n")}
+ list(market.returns.purged=market.returns.purged, innov=innov)
+}
+
+# A function that calls do.one.piece, and works through several
+# different periods to provide the right RHS matrix.
+makeX <- function(market.returns, others,
+ switch.to.innov=rep(TRUE, NCOL(others)),
+ market.returns.purge=TRUE,
+ nlags=5,
+ dates=NULL,
+ verbose=FALSE) {
+ if (verbose) {cat("0. Checking args\n")}
+ stopifnot(all.equal(index(market.returns), index(others)),
+ length(switch.to.innov)==NCOL(others))
+ if (!is.null(dates)) {
+ stopifnot(class(dates) == "Date")
+ }
+ if (verbose) {cat("1. Checking dates.\n")}
+ if (is.null(dates)) {
+ dates <- c(start(market.returns),end(market.returns))
+ }
+ if(head(dates,1)!=head(index(market.returns),1)){
+ stop("Start date provided and the start date of the dataset do not match \n")
+ }
+ if(tail(dates,1)!=tail(index(market.returns),1)){
+ stop("End date provided and the end date of the dataset do not match \n")
+ }
+ if (verbose) {cat("2. Run through all the pieces --\n")}
+ for (i in 1:(length(dates)-1)) {
+ t1 <- dates[i]
+ t2 <- dates[i+1]
+ if (i != (length(dates)-1)) {t2 <- t2 -1}
+ if (verbose) {
+ cat(" Focusing down from date = ", as.character(t1), " to ", as.character(t2), "\n")
+ }
+ tmp.market.returns <- window(market.returns, start=t1, end=t2)
+ tmp.others <- window(others, start=t1, end=t2)
+ a <- do.one.piece(tmp.market.returns, tmp.others, switch.to.innov, market.returns.purge, nlags, verbose)
+ if (i > 1) {
+ res.market.returns <- c(res.market.returns, a$market.returns.purged)
+ res.innov <- rbind(res.innov, a$innov)
+ } else {
+ res.market.returns <- a$market.returns.purged
+ res.innov <- a$innov
+ }
+ }
+ if (verbose) {cat("2. Make a clean X and send it back --\n")}
+ X <- cbind(res.market.returns, res.innov)
+ if (NCOL(res.innov) == 1) {colnames(X) <- c("market.returns","z")}
+ else {colnames(X) <- c("market.returns", colnames(res.innov))}
+ X
+}
+
+
+############################################
+## Summary, print and plot functions for AMM
+############################################
+summary.amm <- function(amm) {
+ cat("\n", "Summary statistics of exposure: \n")
+ sstats <- cbind(amm$exposure, amm$s.exposure,
+ amm$exposure/amm$s.exposure)
+ colnames(sstats) <- c("Exposure", "Std.Err", "t statistic")
+ rownames(sstats) <- names(amm$exposures)
+ print(sstats)
+ cat("\n","Linear model AMM results: ","\n"); class(amm) <- "lm"; print(summary(amm))
+}
+
+print.amm <- function(amm){
+ cat("\n")
+ print(amm$call)
+ cat("\n","Coefficients:","\n")
+ print(amm$coef)
+ cat("\n","Exposures:","\n")
+ print(amm$exposures)
+}
+
+plot.amm <- function(amm){
+ tmp.x <- zoo(as.numeric(resid(amm)), as.Date(names(resid(amm))))
+ tmp.f <- zoo(amm$model$firm.returns, index(tmp.x))
+ tmp <- merge(tmp.x,tmp.f)
+ colnames(tmp) <- c("amm.residuals","firm.returns")
+ plot(tmp, screen=1, lty=1:2, lwd=2, col=c("indian red", "navy blue"),ylab="",
+ xlab="")
+ legend("topleft",legend=c("AMM residual","Firm returns"),lty=1:2, lwd=2,
+ col=c("indian red", "navy blue"), bty='n')
+}
Modified: todo.org
===================================================================
--- todo.org 2014-01-22 10:01:36 UTC (rev 161)
+++ todo.org 2014-02-06 12:50:39 UTC (rev 162)
@@ -4,7 +4,7 @@
- should work with data.frames, zoo, xts
- create a generic "eventstudies" class
- DESIGN!!!
-
+
* phys2eventtime
- more testing: all output is unitmissing
- parallelise (provide an argument to the user)
@@ -12,14 +12,14 @@
> stopifnot(sum(outcomes=="success") == NCOL(z.e))
- if data is missing within the width around the event date, we
simply discard that particular observation from the dataset!
-
+
* AMM()
- vectorise it
- parallelise it
- don't ask the user for amm.type in the eventstudies() function
- separate vignette
- Rewriting? # XXX
-
+
* Vignette
- improvements in AMM()
- missing entries in .bib
More information about the Eventstudies-commits
mailing list