global factor data example

[timelyportfolio]

https://jkpfactors.com/

webr::install("curl", repos = "https://timelyportfolio.github.io/webr_repo/")
webr::install("rolloptim", repos = c("https://timelyportfolio.github.io/webr_repo/", 'https://repo.r-wasm.org', 'https://cran.r-universe.dev'))
webr::install("FactorAnalytics", repos = c("https://timelyportfolio.github.io/webr_repo/", 'https://repo.r-wasm.org', 'https://cran.r-universe.dev'))
library(quantmod)
library(rolloptim)
library(PortfolioAnalytics)
library(FactorAnalytics)
library(DEoptim)
library(GenSA)
library(ROI)
library(ROI.plugin.quadprog)
library(ROI.plugin.deoptim)
temp <- tempfile(fileext=".zip")
download.file("https://jkpfactors.com/data/[usa]_[all_themes]_[monthly]_[vw_cap].zip", temp)
unzip(temp)

factors <- read.csv("/home/web_user/[usa]_[all_themes]_[monthly]_[vw_cap].csv")
factors_xts <- na.omit(do.call(
  merge, 
  lapply(split(factors, factors$name), function(f) {
    x <- xts(
      f$ret,
      order.by = as.Date(f$date)
    )
    colnames(x) <- f$name[[1]]
    x
  })
))

plot.zoo(factors_xts)
chart.CumReturns(factors_xts, legend.loc="topleft")
#chart.RiskReturnScatter(factors_xts)

# rolloptim appears to result in 0 weights in example
n_vars <- 3
n_obs <- 15
x <- matrix(rnorm(n_obs * n_vars), nrow = n_obs, ncol = n_vars)
y <- rnorm(n_obs)
mu <- roll::roll_mean(x, 5)
xx <- roll::roll_crossprod(x, x, 5)
xy <- roll::roll_crossprod(x, y, 5)
sigma <- roll::roll_cov(x, width = 5)
# rolling optimization to minimize variance
roll_min_var(sigma)
# rolling optimization to maximize mean
roll_max_mean(mu)
# rolling optimization to minimize residual sum of squares
roll_min_rss(xx, xy)
# rolling optimization to maximize utility
roll_max_utility(mu, sigma, lambda = 1)

# try PortfolioAnalytics
pspec <- portfolio.spec(assets=colnames(factors_xts))
pspec <- add.constraint(portfolio=pspec,
  type="weight_sum",
  min_sum=0.99,
  max_sum=1.01
) %>%
  add.constraint(type = "box", min = 0.0, max = 0.25) %>%
  add.objective(
    type='risk',
    name='ETL',
    arguments=list(p=0.95)
  ) %>%
  add.objective(
    type='return',
    name='mean'
  )

# Rglpk is not available in webr :(
opt_port_deoptim <- optimize.portfolio(
  na.omit(factors_xts),
  pspec,
  optimize_method = "DEoptim",
  trace = TRUE
)
opt_port_gensa <- optimize.portfolio(
  na.omit(factors_xts),
  pspec,
  optimize_method = "GenSA",
  trace = TRUE
)

#optimize.portfolio(
#  na.omit(factors_xts),
#  pspec,
#  optimize_method = "random"
#)

plot(opt_port_gensa, chart.assets=TRUE, risk.col = "ETL")
chart.RiskReward(opt_port_gensa, return.col="mean", risk.col="ETL",
  chart.assets=TRUE, xlim=c(0.01, 0.05), main="Maximum Return")

portf_minvar <- portfolio.spec(assets=colnames(factors_xts))
# Add full investment constraint to the portfolio object
portf_minvar <- add.constraint(portfolio=portf_minvar, type="full_investment")
# Add objective to minimize variance
portf_minvar <- add.objective(portfolio=portf_minvar, type="risk", name="var")
opt_port_minvar <- optimize.portfolio(
  R = na.omit(factors_xts),
  portfolio = portf_minvar,
  optimize_method = "ROI",
  trace = TRUE
)
plot(opt_port_minvar, chart.assets=TRUE, risk.col = "var")
chart.RiskReward(opt_port_minvar, return.col="mean", risk.col="var",
  chart.assets=TRUE, main="Minimum Variance")

bt_gmv <- optimize.portfolio.rebalancing(
  R = na.omit(factors_xts),
  portfolio = portf_minvar,
  optimize_method = "ROI",
  rebalance_on = "quarters",
  training_period = 36
)

opt_port_meanvar <- optimize.portfolio(
  R = na.omit(factors_xts),
  portfolio = portf_minvar %>%
    add.objective(
      type='return',
      name='mean'
    ) %>%
    add.constraint(
      type="box",
      min=0,
      max=0.4
    ),
  optimize_method = "ROI",
  trace = TRUE
)
plot(opt_port_meanvar, chart.assets=TRUE, risk.col = "var")

weights_meanvar <- as.vector(extractWeights(opt_port_meanvar))
# calculate portfolio returns
port_rets <- apply(
  weights_meanvar * factors_xts,
  1,
  sum
)
port_xts <- xts(port_rets, order.by = index(factors_xts))
colnames(port_xts) <- "meanvar_port"
charts.PerformanceSummary(port_xts)
all_xts <- merge(port_xts,factors_xts)
df <- data.frame(all_xts)
df$date <- index(all_xts)
lm(meanvar_port~.,data=df)
fit <- FactorAnalytics::fitTsfm(asset.names = colnames(all_xts[,1]),
               factor.names = colnames(all_xts[,-1]), 
               data=all_xts)
summary(fit)
fitted(fit)

[timelyportfolio]

library(PerformanceAnalytics)
library(latticeExtra)
webr::install("curl", repos = "https://timelyportfolio.github.io/webr_repo/")
library(quantmod)

temp <- tempfile(fileext=".zip")
download.file("https://jkpfactors.com/data/[usa]_[all_themes]_[monthly]_[vw_cap].zip", temp)
unzip(temp)

factors <- read.csv("/home/web_user/[usa]_[all_themes]_[monthly]_[vw_cap].csv")
factors_xts <- na.omit(do.call(
  merge, 
  lapply(split(factors, factors$name), function(f) {
    x <- xts(
      f$ret,
      order.by = as.Date(f$date)
    )
    colnames(x) <- f$name[[1]]
    x
  })
))

PerformanceAnalytics::chart.Correlation(factors_xts)
# lag.xts causes problems in webr wasm
#   so use xts:::Lag.xts instead
roc_Lag_xts <- function(x, n = 1) {
  x / xts:::Lag.xts(x, n) - 1
}
latticeExtra::horizonplot(roc_Lag_xts(cumprod(1+factors_xts), n = 12), origin = 0, horizonscale = 0.1)

#make a function
#get cumulative of individual components as price
#endpoints will determine periodic rebalancing
#so if you want something other than monthly
calcrebal <- function(data, rebal.period = "months") {
  data.cumul <- data.frame(
    lapply(
      cumprod(1+data),
      FUN=function(x){
        x[endpoints(x,on=rebal.period),]
      }
    )
  )
  data.cumul <- as.xts(
    data.cumul,
    orderBy=as.Date(
      index(
        data[endpoints(data,on=rebal.period),]
      )
    )
  )
  colnames(data.cumul) <- colnames(data)
  #get the returns for a non rebalanced portfolio
  #starting point will be cumulative return for each by themselves
  #so just divide all the monthly values by the beginning value
  bh.cumul <-
    data.cumul /
    matrix(
      data.cumul[1,],
      nrow=NROW(data.cumul),
      ncol=NCOL(data.cumul),
      byrow=TRUE
    )
  #test our calculation graphically
  #should look exactly the same except scale
  #plot.zoo(merge(bh.cumul,data.cumul),nc=2)
  #now let's calculate cumulative at the portfolio level
  #multiple all by 1/N or 1/ncol for equal-weighting
  #then sum by row
  portfolio <- list()
  portfolio$bh <- as.xts(
    apply(
      bh.cumul * 1/NCOL(bh.cumul),
      MARGIN = 1,
      FUN = sum
    ),
    orderBy = as.Date(index(bh.cumul))
  )
  #get the returns for a monthly rebalanced portfolio
  #since we are looking at monthly, we get monthly returns
  #then multiply each by 1/NCOL then sum returns by row
  #get monthly returns
  rebal.cumul <- roc_Lag_xts(data.cumul,1)
  #make first 0 instead of NA to start at 1
  rebal.cumul[1,] <- 0
  portfolio$rebal <- as.xts(
    cumprod(apply(
      rebal.cumul * 1/NCOL(rebal.cumul),
      MARGIN = 1,
      FUN = sum) + 1
    ),
    orderBy = as.Date(index(rebal.cumul))
  )
  #get all indexes in same format,
  #same class, etc. so merge will be proper
  index(portfolio$bh) <- 
    index(portfolio$rebal) <- 
    index(data.cumul)
  return (list(data=data.cumul,portfolio=portfolio))
}
#use our function with our morningstar us style categories
# category_rebal <- calcrebal(ret_xts[,c(3,4,5,6,7,12,13)],"months")
# category_rebal <- calcrebal(ret_xts[,c(5,12)],"months") # low risk and size
category_rebal <- calcrebal(factors_xts[,c(3,13,7)],"months") # value, investment, profit growth
# category_rebal <- calcrebal(ret_xts,"months")

plotdata <- merge(
  log(category_rebal$data),
  log(category_rebal$portfolio$bh),
  log(category_rebal$portfolio$rebal),
  #best component - rebal portfolio
  apply(
    category_rebal$data,
    MARGIN=1,FUN=max
  )-category_rebal$portfolio$rebal,
  #rebalancing bonus(Rb)
  category_rebal$portfolio$rebal-category_rebal$portfolio$bh
)
colnames(plotdata) <- c(
  colnames(category_rebal$data),
  "buyhold",
  "rebal",
  "bestminusrebal",
  "rebalbonus"
)

# xtsExtra::plot.xts(  
#  plotdata,
#   #screens=1,
#   screens=c(rep(1,NCOL(category_rebal$data)),2,2,3,4),
#   layout.screens=matrix(c(1,2,1,2,3,4),ncol=2,byrow=TRUE),
#   ylim = matrix(
#     c(
#       -1,10,
#       -1,10,
#       -2,200,
#       -3,10),ncol=2,byrow=TRUE),
#   auto.legend=TRUE,
#   main = "Rebalancing Bonus (Nardon/Kiskiras)\nKelley US Factor Themes"
# )

# use rebal portfolio for relative performance evaluation
r <- roc_Lag_xts(category_rebal$portfolio$rebal,n=1)
index(r) <- as.Date(index(r))
chart.RelativePerformance(factors_xts,r)
chart.RollingRegression(
  factors_xts,r,width=60,attribute="Alpha",na.pad = FALSE, legend.loc ="top", colorset = RColorBrewer::brewer.pal(name="Set2",n=8),
  main = "US Factor Themes Rolling 60-month Alpha to Equal Weight"
)

# momentum strategy buy for width best performing over width style
#ret_xts <- ret_xts["1970::",]
width <- 12
roc <- na.omit(
  rollapplyr(data = factors_xts, FUN = function(x){tail(cumprod(1+x)-1,1)}, width=width, by=width)
)
alpha <- na.omit(
  rollapplyr(data = factors_xts, FUN = function(x){CAPM.jensenAlpha(x,r[index(x),])}, width=width, by=width)
)
maxstyle <- unname(unlist(apply(roc,FUN=which.max, MARGIN=1)))
minstyle <- unname(unlist(apply(roc,FUN=which.min, MARGIN=1)))
alphamax <- unname(unlist(apply(alpha,FUN=which.max, MARGIN=1)))
alphamin <- unname(unlist(apply(alpha,FUN=which.min, MARGIN=1)))

sys_max <- xts(
  mapply(
    function(r,x) {
      r[[1]][x]
    },
    r = apply(X=roc[-1,], MARGIN=1, function(x){list(unname(unlist(x)))}),
    x = maxstyle[-length(maxstyle)],
    USE.NAMES = FALSE
  ),
  order.by = index(roc[-1,])
)
sys_min <- xts(
  mapply(
    function(r,x) {
      r[[1]][x]
    },
    r = apply(X=roc[-1,], MARGIN=1, function(x){list(unname(unlist(x)))}),
    x = minstyle[-length(minstyle)],
    USE.NAMES = FALSE
  ),
  order.by = index(roc[-1,])
)
sys_alpha_max <- xts(
  mapply(
    function(r,x) {
      r[[1]][x]
    },
    r = apply(X=roc[-(1:2),], MARGIN=1, function(x){list(unname(unlist(x)))}),
    x = alphamax[-length(alphamax)],
    USE.NAMES = FALSE
  ),
  order.by = index(roc[-(1:2),])
)
sys_alpha_min <- xts(
  mapply(
    function(r,x) {
      r[[1]][x]
    },
    r = apply(X=roc[-(1:2),], MARGIN=1, function(x){list(unname(unlist(x)))}),
    x = alphamin[-length(alphamin)],
    USE.NAMES = FALSE
  ),
  order.by = index(roc[-(1:2),])
)
charts.PerformanceSummary(
  na.omit(
    merge(
      roc_Lag_xts(cumprod(1+factors_xts),n=width),
      sys_max,
      sys_min,
      sys_alpha_max,
      sys_alpha_min
    )
  ),
  colorset = RColorBrewer::brewer.pal(name="Set3",n=12)
)
chart.CumReturns(
  na.omit(
    merge(
      roc_Lag_xts(cumprod(1+factors_xts),n=width),
      sys_max,
      sys_min,
      sys_alpha_max,
      sys_alpha_min#,
      # rebal portfolio with all
      #roc_Lag_xts(cumprod(1+na.omit(r)),n=width)
    )
  ),
  legend.loc = "topleft",
  colorset = RColorBrewer::brewer.pal(name="Set3",n=12)
)
#chart.RiskReturnScatter(na.omit(merge(roc_Lag_xts(cumprod(1+ret_xts),n=width), sys)))

# monthly system (by=1)
width <- 10
roc2 <- na.omit(
  rollapplyr(data = factors_xts, FUN = function(x){tail(cumprod(1+x)-1,1)}, width=width, by=1)
)
alpha2 <- na.omit(
  rollapplyr(data = factors_xts, FUN = function(x){CAPM.jensenAlpha(x,r[index(x),])}, width=width, by=1)
)
maxstyle2 <- unname(unlist(apply(roc2,FUN=which.max, MARGIN=1)))
minstyle2 <- unname(unlist(apply(roc2,FUN=which.min, MARGIN=1)))
alphamax2 <- unname(unlist(apply(alpha2,FUN=which.max, MARGIN=1)))
alphamin2 <- unname(unlist(apply(alpha2,FUN=which.min, MARGIN=1)))
sys_max2 <- xts(
  mapply(
    function(r,x) {
      r[[1]][x]
    },
    r = apply(X=factors_xts[paste0(range(index(roc2)),collapse="::"),][-1,], MARGIN=1, function(x){list(unname(unlist(x)))}),
    x = maxstyle2[-length(maxstyle2)],
    USE.NAMES = FALSE
  ),
  order.by = index(factors_xts[paste0(range(index(roc2)),collapse="::"),][-1,])
)
sys_min2 <- xts(
  mapply(
    function(r,x) {
      r[[1]][x]
    },
    r = apply(X=factors_xts[paste0(range(index(roc2)),collapse="::"),][-1,], MARGIN=1, function(x){list(unname(unlist(x)))}),
    x = minstyle2[-length(minstyle2)],
    USE.NAMES = FALSE
  ),
  order.by = index(factors_xts[paste0(range(index(roc2)),collapse="::"),][-1,])
)
sys_alpha_max2 <- xts(
  mapply(
    function(r,x) {
      r[[1]][x]
    },
    r = apply(X=factors_xts[paste0(range(index(alpha2)),collapse="::"),][-1,], MARGIN=1, function(x){list(unname(unlist(x)))}),
    x = alphamax2[-length(alphamax2)],
    USE.NAMES = FALSE
  ),
  order.by = index(factors_xts[paste0(range(index(alpha2)),collapse="::"),][-1,])
)
sys_alpha_min2 <- xts(
  mapply(
    function(r,x) {
      r[[1]][x]
    },
    r = apply(X=factors_xts[paste0(range(index(alpha2)),collapse="::"),][-1,], MARGIN=1, function(x){list(unname(unlist(x)))}),
    x = alphamin2[-length(alphamin2)],
    USE.NAMES = FALSE
  ),
  order.by = index(factors_xts[paste0(range(index(alpha2)),collapse="::"),][-1,])
)
chart.CumReturns(
  na.omit(
    merge(
      # roc_Lag_xts(cumprod(1+ret_xts),n=width),
      factors_xts,
      sys_max2,
      sys_min2,
      sys_alpha_max2,
      sys_alpha_min2#,
      # rebal portfolio with all
      #roc_Lag_xts(cumprod(1+na.omit(r)),n=width)
    )
  ),
  legend.loc = "topleft",
  colorset = RColorBrewer::brewer.pal(name="Set3",n=12)
)

[timelyportfolio]

cleaner version of first example

webr::install("curl", repos = "https://timelyportfolio.github.io/webr_repo/")
library(quantmod)
library(PortfolioAnalytics)
library(ROI)
library(ROI.plugin.quadprog)

# use factor data from https://jkpfactors.com/
# Jensen, T., Kelly, B., and Pedersen, L. “Is There a Replication Crisis in Finance?” Journal of Finance (2023)
temp <- tempfile(fileext=".zip")
download.file("https://jkpfactors.com/data/[usa]_[all_themes]_[monthly]_[vw_cap].zip", temp)
unzip(temp)
factors <- read.csv("/home/web_user/[usa]_[all_themes]_[monthly]_[vw_cap].csv")
factors_xts <- na.omit(do.call(
  merge, 
  lapply(split(factors, factors$name), function(f) {
    x <- xts(
      f$ret,
      order.by = as.Date(f$date)
    )
    colnames(x) <- f$name[[1]]
    x
  })
))

plot.zoo(factors_xts)
# chart.CumReturns(factors_xts, legend.loc="topleft")

# try PortfolioAnalytics
# min variance
opt_port_minvar <- optimize.portfolio(
  R = na.omit(factors_xts),
  portfolio = portfolio.spec(assets=colnames(factors_xts)) |>
    add.constraint(type="full_investment") |>
    add.objective(type="risk", name="var"),
  optimize_method = "ROI",
  trace = TRUE
)
plot(opt_port_minvar, chart.assets=TRUE, risk.col = "var")

bt_gmv <- optimize.portfolio.rebalancing(
  R = na.omit(factors_xts),
  portfolio = portf_minvar,
  optimize_method = "ROI",
  rebalance_on = "quarters",
  training_period = 36
)

# mean variance portfolio
opt_port_meanvar <- optimize.portfolio(
  R = na.omit(factors_xts),
  portfolio = portfolio.spec(assets=colnames(factors_xts)) |>
    add.constraint(type="full_investment") |>
    add.objective(type="risk", name="var") |>
    add.objective(type='return', name='mean') |>
    add.constraint(type="box", min=0, max=0.4),
  optimize_method = "ROI",
  trace = TRUE
)
plot(opt_port_meanvar, chart.assets=TRUE, risk.col = "var")

# calculate portfolio returns for mean variance
weights_meanvar <- as.vector(extractWeights(opt_port_meanvar))
port_rets <- apply(
  weights_meanvar * factors_xts,
  1,
  sum
)
port_xts <- xts(port_rets, order.by = index(factors_xts))
colnames(port_xts) <- "meanvar_port"
charts.PerformanceSummary(port_xts)