[Eventstudies-commits] r277 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 2 18:25:49 CEST 2014
Author: chiraganand
Date: 2014-04-02 18:25:48 +0200 (Wed, 02 Apr 2014)
New Revision: 277
Added:
pkg/R/lmAMM.R
Removed:
pkg/R/lmAmm.R
Log:
Renamed lmAMM R file, making consistent with the function name.
Copied: pkg/R/lmAMM.R (from rev 275, pkg/R/lmAmm.R)
===================================================================
--- pkg/R/lmAMM.R (rev 0)
+++ pkg/R/lmAMM.R 2014-04-02 16:25:48 UTC (rev 277)
@@ -0,0 +1,362 @@
+###############
+## One firm AMM
+###############
+## This function takes care of the structural break dates introduced by calculating exposure for sub periods differently
+## This function is used when date are provided in the function
+subperiod.lmAMM <- function(firm.returns,X,nlags=1,verbose=FALSE,dates=NULL,residual=TRUE){
+ ## Creating empty frames
+ if(is.null(dates)){
+ dates.no <- c(start(firm.returns),end(firm.returns))
+ } else{
+ dates.no <- dates
+ }
+ exposures <- data.frame(matrix(NA,ncol=ncol(X),nrow=(length(dates.no)-1)))
+ colnames(exposures) <- colnames(X)
+ sds <- exposures
+ periodnames <- NULL
+
+ ## Getting firm exposure, amm residuals
+ if(is.null(dates)){
+ res <- lmAMM(firm.returns,X,verbose=verbose,nlags=nlags)
+ if(is.null(res)!=TRUE){
+ exposures <- res$exposure
+ sds <- res$s.exposure
+ m.residuals <- xts(res$residuals,as.Date(attr(res$residuals,"names")))
+ if(residual==TRUE){
+ m.residuals <- xts(res$residuals,as.Date(attr(res$residuals,"names")))
+ }
+ rval <- list(exposures=exposures,sds=sds,residuals=m.residuals)
+ } else {
+ rval <- NULL
+ }
+ }else{
+ tmp <- window(firm.returns,start=dates[1],end=dates[1+1])
+ rhs <- window(X,start=dates[1],end=dates[1+1])
+ res <- lmAMM(firm.returns=tmp,
+ X=rhs,
+ verbose=verbose,
+ nlags=nlags)
+ exposures[1,] <- res$exposure
+ periodnames <- c(periodnames,paste(dates[1],dates[1+1],sep=" TO "))
+ sds[1,] <- res$s.exposure
+ m.residuals <- xts(res$residuals,as.Date(attr(res$residuals,"names")))
+ colnames(m.residuals) <- paste(dates[1],"to",dates[1+1],sep=".")
+ for(i in 2:(length(dates)-1)){
+ tmp <- window(firm.returns,start=dates[i],end=dates[i+1])
+ rhs <- window(X,start=dates[i],end=dates[i+1])
+ res <- lmAMM(firm.returns=tmp,
+ X=rhs,
+ verbose=verbose,
+ nlags=nlags)
+ exposures[i,] <- res$exposure
+ periodnames <- c(periodnames,paste(dates[i],dates[i+1],sep=" TO "))
+ sds[i,] <- res$s.exposure
+ period.resid <- xts(res$residuals,as.Date(attr(res$residuals,"names")))
+ colnames(period.resid) <- paste(dates[i],"to",dates[i+1],sep=".")
+ m.residuals <- merge(m.residuals, period.resid, all=TRUE)
+ }
+ rownames(exposures) <- rownames(sds) <- periodnames
+ rval <- list(exposures=exposures,sds=sds,residuals=m.residuals)
+ }
+ return(rval)
+}
+
+########################
+# Many firms AMM
+########################
+manyfirmssubperiod.lmAMM <-
+function(firm.returns,X,
+ lags,dates=NULL, periodnames=NULL,verbose=FALSE){
+ if(is.null(dates)){
+ dates=c(start(X),end(X))
+ 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(firm.returns)
+
+ # Let's get "exposure' and 'sds'. Setting up structures:-
+
+ exposures <- matrix(NA,nrow=nfirms,ncol=nperiods*ncol(X))
+ exposures <- as.data.frame(exposures)
+ rownames(exposures) <- colnames(firm.returns)
+ tmp <- NULL
+ for(i in 1:length(periodnames)){
+ for(j in 1:NCOL(X)){
+ tmp <- c(tmp, paste(colnames(X)[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(X)),
+ s.exposures=rep(NA,ncol(X)))
+
+ for(i in 1:NCOL(firm.returns)){
+ cat("AMM estimation for",colnames(firm.returns)[i],"\n")
+ if (verbose) {cat ("AMM estimation for", colnames(firm.returns)[i], "\n")}
+ stock.return <- firm.returns[,i]
+ dataset <- cbind(stock.return, X) # 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(firm.returns)[i],] <- this.exp
+ sds[colnames(firm.returns)[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)
+}
+
+
+###########################
+# Maintaining 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(object, ...) {
+ cat("\n", "Summary statistics of exposure: \n")
+ sstats <- cbind(object$exposure, object$s.exposure,
+ object$exposure/object$s.exposure)
+ colnames(sstats) <- c("Exposure", "Std.Err", "t statistic")
+ rownames(sstats) <- names(object$exposures)
+ print(sstats)
+ cat("\n","Linear model AMM results: ","\n");
+ class(object) <- "lm";
+ print.default(summary.default(object))
+}
+
+print.amm <- function(x, ...){
+ cat("\n")
+ print(x$call)
+ cat("\n","Coefficients:","\n")
+ print(x$coef)
+ cat("\n","Exposures:","\n")
+ print.default(x$exposures)
+}
+
+plot.amm <- function(x, ...){
+ tmp.x <- zoo(as.numeric(resid(x)), as.Date(names(resid(x))))
+ tmp.f <- zoo(x$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')
+}
Deleted: pkg/R/lmAmm.R
===================================================================
--- pkg/R/lmAmm.R 2014-04-02 16:14:17 UTC (rev 276)
+++ pkg/R/lmAmm.R 2014-04-02 16:25:48 UTC (rev 277)
@@ -1,362 +0,0 @@
-###############
-## One firm AMM
-###############
-## This function takes care of the structural break dates introduced by calculating exposure for sub periods differently
-## This function is used when date are provided in the function
-subperiod.lmAMM <- function(firm.returns,X,nlags=1,verbose=FALSE,dates=NULL,residual=TRUE){
- ## Creating empty frames
- if(is.null(dates)){
- dates.no <- c(start(firm.returns),end(firm.returns))
- } else{
- dates.no <- dates
- }
- exposures <- data.frame(matrix(NA,ncol=ncol(X),nrow=(length(dates.no)-1)))
- colnames(exposures) <- colnames(X)
- sds <- exposures
- periodnames <- NULL
-
- ## Getting firm exposure, amm residuals
- if(is.null(dates)){
- res <- lmAMM(firm.returns,X,verbose=verbose,nlags=nlags)
- if(is.null(res)!=TRUE){
- exposures <- res$exposure
- sds <- res$s.exposure
- m.residuals <- xts(res$residuals,as.Date(attr(res$residuals,"names")))
- if(residual==TRUE){
- m.residuals <- xts(res$residuals,as.Date(attr(res$residuals,"names")))
- }
- rval <- list(exposures=exposures,sds=sds,residuals=m.residuals)
- } else {
- rval <- NULL
- }
- }else{
- tmp <- window(firm.returns,start=dates[1],end=dates[1+1])
- rhs <- window(X,start=dates[1],end=dates[1+1])
- res <- lmAMM(firm.returns=tmp,
- X=rhs,
- verbose=verbose,
- nlags=nlags)
- exposures[1,] <- res$exposure
- periodnames <- c(periodnames,paste(dates[1],dates[1+1],sep=" TO "))
- sds[1,] <- res$s.exposure
- m.residuals <- xts(res$residuals,as.Date(attr(res$residuals,"names")))
- colnames(m.residuals) <- paste(dates[1],"to",dates[1+1],sep=".")
- for(i in 2:(length(dates)-1)){
- tmp <- window(firm.returns,start=dates[i],end=dates[i+1])
- rhs <- window(X,start=dates[i],end=dates[i+1])
- res <- lmAMM(firm.returns=tmp,
- X=rhs,
- verbose=verbose,
- nlags=nlags)
- exposures[i,] <- res$exposure
- periodnames <- c(periodnames,paste(dates[i],dates[i+1],sep=" TO "))
- sds[i,] <- res$s.exposure
- period.resid <- xts(res$residuals,as.Date(attr(res$residuals,"names")))
- colnames(period.resid) <- paste(dates[i],"to",dates[i+1],sep=".")
- m.residuals <- merge(m.residuals, period.resid, all=TRUE)
- }
- rownames(exposures) <- rownames(sds) <- periodnames
- rval <- list(exposures=exposures,sds=sds,residuals=m.residuals)
- }
- return(rval)
-}
-
-########################
-# Many firms AMM
-########################
-manyfirmssubperiod.lmAMM <-
-function(firm.returns,X,
- lags,dates=NULL, periodnames=NULL,verbose=FALSE){
- if(is.null(dates)){
- dates=c(start(X),end(X))
- 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(firm.returns)
-
- # Let's get "exposure' and 'sds'. Setting up structures:-
-
- exposures <- matrix(NA,nrow=nfirms,ncol=nperiods*ncol(X))
- exposures <- as.data.frame(exposures)
- rownames(exposures) <- colnames(firm.returns)
- tmp <- NULL
- for(i in 1:length(periodnames)){
- for(j in 1:NCOL(X)){
- tmp <- c(tmp, paste(colnames(X)[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(X)),
- s.exposures=rep(NA,ncol(X)))
-
- for(i in 1:NCOL(firm.returns)){
- cat("AMM estimation for",colnames(firm.returns)[i],"\n")
- if (verbose) {cat ("AMM estimation for", colnames(firm.returns)[i], "\n")}
- stock.return <- firm.returns[,i]
- dataset <- cbind(stock.return, X) # 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(firm.returns)[i],] <- this.exp
- sds[colnames(firm.returns)[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)
-}
-
-
-###########################
-# Maintaining 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(object, ...) {
- cat("\n", "Summary statistics of exposure: \n")
- sstats <- cbind(object$exposure, object$s.exposure,
- object$exposure/object$s.exposure)
- colnames(sstats) <- c("Exposure", "Std.Err", "t statistic")
- rownames(sstats) <- names(object$exposures)
- print(sstats)
- cat("\n","Linear model AMM results: ","\n");
- class(object) <- "lm";
- print.default(summary.default(object))
-}
-
-print.amm <- function(x, ...){
- cat("\n")
- print(x$call)
- cat("\n","Coefficients:","\n")
- print(x$coef)
- cat("\n","Exposures:","\n")
- print.default(x$exposures)
-}
-
-plot.amm <- function(x, ...){
- tmp.x <- zoo(as.numeric(resid(x)), as.Date(names(resid(x))))
- tmp.f <- zoo(x$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')
-}
More information about the Eventstudies-commits
mailing list