Time-series-momentum.Rmd

---
title: 'Time Series Momentum, a replication'
subtitle: "(draft)"
author: "Vito Lestingi, Justin M Shea"
date: "`r Sys.Date()`"
abstract: "This document includes replication material on some academic and practitioners' literature instrumental for our RGSoC 2020 project. The document itself is meant to be completely reproducible."
output:
  rmarkdown::pdf_document:
    citation_package: biblatex
bibliography: references.bib
vignette: >
  %\VignetteIndexEntry{"Time Series Momentum"}
  %\VignetteEngine{rmarkdown::render}
  %\VignetteEncoding{UTF-8}
---


```{r, include=FALSE}
options(tinytex.verbose = TRUE)
```


```{r setup, include = FALSE}
knitr::opts_knit$set(
  root.dir = normalizePath('..')
)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo = TRUE,
  warning = FALSE,
  error = FALSE
)

# Load libraries
library(xts)
library(plm)
library(FactorAnalytics)
```

```{r, echo=FALSE}
# Mentors, please feel free to add yourself to the authors' field if you wish.
```

# Introduction
\textcite{moskowitz-ooi-pedersen-2012} (MOP hereafter) study an asset # anomaly they name *time series momentum*, which is related but different from the *momentum* effect in that the latter has a cross-sectional relative nature with respect to assets clusters while the former is directly linked with single assets returns. They find the anomaly consistent both across different asset classes and markets. Also, they confirm this effect to be robust among more illiquid instruments.

# Data and methodology
TODO: Present data series we used, compare with analogies and differences with respect to data authors used.

# The cross-section of time series momentum
The excess returns 
$$
\frac{r_t^s}{\sigma_{t-1}^s} = \alpha + \beta_h\frac{r_{t-h}^s}{\sigma_{t-h-1}^s} + \epsilon_t^s
$$

where assets returns are scaled by their ex-ante volatility $\sigma_{t-1}^s$, with annualized variance being
$$
\sigma_t^2 = 261\sum_{i=0}^{\infty}(1 - \delta)\delta^i(r_{t-i-1} - \bar{r}_t)^2
$$
$\bar{r}_t$ returns exponentially weighted moving average (EWMA) and weight $\delta$ so that $\delta/(1 - \delta) = 60$ days.
Returns of TSMOM trading strategies on the instrument $s$ at month $t$ are systematically obtained considering the excess returns sign of $s$ over the past $k$ months and then acquiring or selling the instrument during the subsequent $h$ months. These two parameters determine a family of TSMOM strategies and are called *lookback period* and *holding period*, respectively.  
It follows that the strategy return at time $t$, defined $r_t^{\textrm{TSMOM}(k,h)}$, represents the average return across all instrument portfolios at that time, i.e. the return on the portfolio that was constructed in all observable past months. These returns are then average across all instruments include or within each asset class.

To explain these returns and assess whether they held abnormal performance, \textcite{moskowitz-ooi-pedersen-2012} study the regression specification
$$
r_t^{\textrm{TSMOM}(k,h)} = \alpha + \beta_{1}MSCI_{t} + \beta_{2}GSCI_{t} + \beta_{3}BOND_{t} + \beta_{4}SMB_{t} + \beta_{5}HML_{t} + \beta_{6}UMD_{t} + \epsilon_{t}
$$
where $MSCI$ is the MSCI World Index, $GSCI$ the S&P Goldman Sachs Commodity Index, $BOND$ is Barclay's Aggregate Bond Index ([Bloomberg Barclays Global Aggregate Index](https://data.bloomberglp.com/indices/sites/2/2016/08/Factsheet-Global-Aggregate.pdf) at the time of writing) and $SMB$, $HML$, $UMD$ are the usual Fama-French-Carhart factors. In what follows, given data series availability from authors, we study the TSMOM strategy with a 12 months lookback period and holding period of a month, that is $r_t^{\textrm{TSMOM}(12,1)}$, for each asset class and in the all assets aggregate.

```{r Load Factors and Portfolios Returns Data}
# TSM portfolios
path.parser <- system.file('parsers', 
                           'TSM.R', 
                           package = 'ExpectedReturns')
source(path.parser)
TSM <- xts::xts(TSM[, -1], order.by=TSM$DATE)
# MSCI
# NOTE: returns in decimal unit
path.parser <- system.file('parsers', 
                           'MSCI-WI.R', 
                           package = 'ExpectedReturns')
source(path.parser)
# NOTE: 
# - in place of 'GSCI' we use 'CM.MARKET'
# - in place of 'BOND' we use 'FI.MARKET'
path.parser <- system.file('parsers', 
                           'CFP.R', 
                           package = 'ExpectedReturns')
source(path.parser)
CM.MARKET <- CFP[, 'CM.MARKET']
FI.MARKET <- CFP[, 'FI.MARKET']
# VME
path.parser <- system.file('parsers', 
                           'VME-Factors.R', 
                           package = 'ExpectedReturns')
source(path.parser)
VME.FACTORS <- VME.Factors[, c('VAL.EVR', 'MOM.EVR')]
# FFC factors

path.parser <- system.file('parsers', 
                           'FFdownloads_factors_package.R', 
                           package = 'ExpectedReturns')
source(path.parser)
## Function to parse through 'ffdownloads' package available monthly data from "http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html" - Bryan
FFfactors <- FFfactors_xts(x)
## assigning same name as used in this and other scripts.
FF3 <- FFfactors[["x_F-F_Research_Data_Factors"]]
MOM <- FFfactors[["x_F-F_Momentum_Factor"]]
#FF3 <- ExpectedReturns::GetFactors('FF3', freq='monthly')
#MOM <- ExpectedReturns::GetFactors('MOM', freq='monthly')
## aded as.Date() wrap since date format in new data comes in as "%b %Y" (May 2022). Need to properly adjust FFfactors function so this is taken care of before.
min.tp <- max(first(as.Date(index(FF3))), first(as.Date(index(MOM))))
max.tp <- min(last(as.Date(index(FF3))), last(as.Date(index(FF3))))
# min.tp <- max(first(index(FF3)), first(index(MOM)))
# max.tp <- min(last(index(FF3)), last(index(MOM)))
days.diff <- diff(seq.Date(min.tp, max.tp, by='month'))[-1]
ff.dates <- c(min.tp, min.tp + cumsum(as.numeric(days.diff)))
# VIX Index

path.parser <- system.file('parsers', 
                           'VIX-FRED.R', 
                           package = 'ExpectedReturns')
source(path.parser)
VIX.RET <- VIX.cls.monthly$VIX.RET
# VIX Index top ~20% extremes
vix20idxs <- order(
  abs(VIX.RET$VIX.RET), decreasing=TRUE
)[1:round(nrow(VIX.RET) * 0.2)]
VIX20 <- VIX.RET[vix20idxs, ]
colnames(VIX20) <- 'VIX.TOP.20'
# TED Spread
path.parser <- system.file('parsers', 
                           'TED-Spread.R', 
                           package = 'ExpectedReturns')
source(path.parser)
colnames(TED.SPREAD) <- 'TED'
# TED Spread top ~20% extremes
ted20idxs <- order(
  abs(TED.SPREAD$TED), decreasing=TRUE
)[1:round(nrow(TED.SPREAD) * 0.2)]
TED20 <- TED.SPREAD[ted20idxs, ]
colnames(TED20) <- 'TED.TOP.20'
```

```{r Prepare data set}
# FFC4 factors
FFC4 <- merge(FF3[ff.dates, ], MOM[ff.dates, ])
# MSCI World Index
data <- merge(MSCI.WI$RET, FFC4)
data$MSCI.RET <- zoo::na.locf(data$RET)
data$MSCI.RF <- data$MSCI.RET - data$RF
data$RET <- NULL

# Bonds factors data
# data <- merge(CRP, data)
# data$GOVT.XS <- na.fill(data$GOVT.XS, c(NA, 'extend', NA))
# data$CORP.XS <- na.fill(data$CORP.XS, c(NA, 'extend', NA))
# Convert data set
# data <- data.frame(
#   DATE=ff.dates,
#   data[ff.dates, c(colnames(CRP), 'MSCI.RF', colnames(FFC4))],
#   row.names=NULL
# )
# tp <- 1:max(which(!is.na(data$CORP.XS)))
# data <- data[tp, ]

# Bonds factor
data <- merge(data, FI.MARKET)
# Commodities Market
data <- merge(data, CM.MARKET)
# VME factors
data <- merge(data, VME.FACTORS)
# VIX Index
data <- merge(data, VIX.RET)
# TED Spread
data <- merge(data, TED.SPREAD)
# TSMOM(12, 1) strategy returns
data <- merge(data, TSM)

# NOTE:
# All series are considered relative to FF dates, which are at month-end.
# Usually a dates mismatch of one day can exist around the month-end, for reasons 
# among which publication date discrepancies or subsequent corrections. 
# When this happens we simply consider last available values with respect to the 
# month-end, as those are dates most series we work with refer to.
data <- zoo::na.locf(data)
# NOTE:
# TED-Spread and VIX top ~20% series are used as they are constructed
data <- merge(data, TED20, VIX20)

data <- data.frame(
  DATE=ff.dates,
  data[ff.dates, ],
  row.names=NULL
)
```

```{r TSM_Table_2_Time-series_Regressions}
# Indexes
date.id <- matrix(1:nrow(data), dimnames=list(NULL, 'DATE.ID'))
data <- cbind(data, date.id)
# The period used in the paper for running the time series regression is Jan 1985 - Dec 2009
# The AQR website TSMOM portfolios from which we sourced our data for this replication start # reporting data from that date, but keep the data up to date. For the purposes of
# replicating the time period in the paper, we will need to truncate the 'data' object after # Dec 2009. The resulting t-stats of the alpha intercepts are more inline with the paper if
# we do. For the updated alpha t-stats the user can comment this line of code. 
data <- data[-((which(data$DATE=="2009-12-31")+1):nrow(data)),]
# Regressions variables
y <- colnames(TSM)
X <- c('MSCI.RF', 'CM.MARKET', 'FI.MARKET', 'SMB', 'HML', 'MOM')
# Time-series regressions
tsmom.ts.reg <- lapply(1:length(y), function(x) {
  model.formula <- formula(
    paste(
      y[x], paste(X, collapse='+'), 
      sep='~'
    )
  )
  plm::plm(
    model.formula, data=data,
    model='pooling', index='DATE.ID'
  )
})
lapply(tsmom.ts.reg, summary)
```

```{r TSM_Table_2_Time-series_Regressions_FA}
#refactored
data.test = data
row.names(data.test) <-data.test$DATE
data.test$DATE <- NULL
tsmom.ts.reg.test <- lapply(1:length(y), function(x) {
  model.formula <- formula(
    paste(
      y[x], paste(X, collapse='+'), 
      sep='~'
    )
  )
  fitTsfm(factor.names=X, asset.names=y[x], data = data.test)
}) 
lapply(tsmom.ts.reg.test, summary)
```

```{r TSM_Table_3_Time-series_Regressions}
# NOTE:
# for the last three models we run regressions on monthly series as opposed to 
# quarterly data.
y <- 'TSMOM' # diversified TSMOM(12,1)
X <- list(
  ffc4.msci=c('MSCI.RF', 'SMB', 'HML', 'MOM')
  , amp3=c('MSCI.RF', 'VAL.EVR', 'MOM.EVR')
  , msci=c('MSCI.RF', 'I(MSCI.RF^2)')
  , ted=c('TED')
  , ted20=c('TED.TOP.20')
  , vix=c('VIX.RET')
  , vix20=c('VIX.TOP.20')
)
tsmom.div.ts.reg <- lapply(X, function(x) {
  model.formula <- formula(
    paste(
      y, paste(x, collapse='+'), 
      sep='~'
    )
  )  
  plm::plm(
    model.formula, data=data,
    model='pooling', index='DATE.ID'
  )
})
lapply(tsmom.div.ts.reg, summary)
```

```{r TSM_Table_3_Time-series_Regressions_FA}
#refactored by Jiarui
data.test = data
row.names(data.test) <-data.test$DATE
data.test$DATE <- NULL
data.test['MSCI.RF.square']<-data.test$MSCI.RF^2 
X <- list(    
  ffc4.msci=c('MSCI.RF', 'SMB', 'HML', 'MOM')
  , amp3=c('MSCI.RF', 'VAL.EVR', 'MOM.EVR')
  , msci=c('MSCI.RF', 'MSCI.RF.square')
  , ted=c('TED')
  , ted20=c('TED.TOP.20')
  , vix=c('VIX.RET')
  , vix20=c('VIX.TOP.20')
) 
tsmom.div.ts.reg.test <- lapply(X, function(x) { 
  fitTsfm(factor.names=x, asset.names='TSMOM', data = data.test) 
}) 
lapply(tsmom.div.ts.reg.test, summary)
```

# Time series momentum factor
Let us consider the TSMOM(12, 1) strategy, aggregating returns across all asset classes holds a portfolio called *diversified TSMOM factor* and expressed as
$$
r_{t,t+1}^{\textrm{TSMOM}} = \frac{1}{S_t}\sum_{s=1}^{S_t}\textrm{sign}(r_{t-12,t}^{s})\frac{40\%}{\sigma_t^s}r_{t,t+1}^{s}
$$
with $S_t$ securities investable at time $t$ and a 40% constant annual volatility is chosen by authors because "it is similar to the risk of an average individual stock" and to "make it easier to intuitively compare our portfolios to other in the literature" as, it is consistent with other factors' volatility once averaged over securities.

# Time series momentum vs. cross-sectional momentum
Follow authors, in this section we compare time series momentum and the cross-sectional momentum of \textcite{asness-moskowitz-pedersen-2013}. Over the comparable sample period, our results are close to the ones authors published. They differ in that often our estimates exhibit heavier loading in direct momentum time-series regressions. In particular, an empirical interpretation of small magnitude signs shifts may be that during the last two to three years sample period analyzed by authors financial markets were rather turbulent. A practical reason for differences is that data series corrections have occurred in authors' updated data sets we are working with. Notwithstanding, t-statistics and $R^2$ generally appear to be in line with authors' results and this may be confirmatory of the stable relation between the two momentum strategies both across asset classes and over time.

```{r TSM_Table_5_Panel_A}
# Time series and Cross-sectional Momentum data
XSMOM <- VME.Factors[, c('MOM.AA', 'MOMLS.VME.COM', 'MOMLS.VME.EQ', 'MOMLS.VME.FI', 'MOMLS.VME.FX', 'MOMLS.VME.US90')]
colnames(XSMOM) <- c('XSMOM.ALL', 'XSMOM.COM', 'XSMOM.EQ', 'XSMOM.FI', 'XSMOM.FX', 'XSMOM.US') # naming consistency
#XSMOM <- xts::xts(XSMOM[, -1], order.by=XSMOM$DATE)
mom.data <- merge(TSM, XSMOM)
mom.data <- zoo::na.locf(mom.data)
mom.data <- mom.data[xts::endpoints(mom.data), ]
mom.data$DATE.ID <- 1:nrow(mom.data)
# To check out on paper sample period
# mom.data <- mom.data['1985/2009', ]
mom.data <- data.frame(
  DATE=index(mom.data),
  mom.data,
  row.names=NULL
)
# Time-series regressions
Y <- rep(colnames(TSM), 2)
X.aac <- colnames(XSMOM)[-1]
X <- c(rep(list(X.aac), 5), colnames(XSMOM)[1], as.list(X.aac[-length(X.aac)]))
mom.tsxs.reg <- lapply(1:length(Y), function(x) {
  model.formula <- formula(
    paste(
      Y[x], paste(X[[x]], collapse='+'), 
      sep='~'
    )
  )
  plm::plm(
    model.formula, data=mom.data,
    model='pooling', index='DATE.ID'
  )
})
mom.tsxs.reg <- lapply(mom.tsxs.reg, summary)
names(mom.tsxs.reg) <- Y
mom.tsxs.reg
```

```{r TSM_Table_5_Panel_A_FA}
#refactored by Jiarui
mom.data.test = mom.data
row.names(mom.data.test) <-mom.data.test$DATE 
mom.data.test$DATE <- NULL
mom.tsxs.reg.test <-lapply(1:length(Y), function(x) {
  fitTsfm(factor.names=X[[x]], asset.names=Y[x], data = mom.data.test) 
}) 
lapply(mom.tsxs.reg.test, summary)
 
mom.tsxs.reg.test <- lapply(mom.tsxs.reg.test, summary)
names(mom.tsxs.reg.test) <- Y
mom.tsxs.reg.test
```

Next, we try to gain insights on what factors better explain time-series momentum.
```{r Table_5_Panel_C}
## Panel C
# NOTE: 'DJCS MF' and 'DJCS MACRO' data missing
data <- merge(
  data, 
  data.frame(
    DATE=XSMOM[, 1], 
    XSMOM[, -1],
    row.names=NULL
  )
)
Y <- c('XSMOM.ALL', 'XSMOM.COM', 'XSMOM.EQ', 'XSMOM.FI', 'XSMOM.FX', 'SMB', 'HML', 'MOM')
tsmom.div.all <- lapply(Y, function(y) {
  plm::plm(
    formula(paste(y, 'TSMOM', sep='~')), 
    data=data, model='pooling', index='DATE.ID'
  )
})
lapply(tsmom.div.all, summary)
```

# References