Optimizing Your Personal Investment Portfolio with ETFs

#download the list of ETFs at: https://www.nasdaq.com/etfs/list

from alpha_vantage.timeseries import TimeSeries
import numpy as np
import pandas as pd
from pandas import DataFrame as df
import os
import csv
import datatable as dt

# Specify the API key
ts = TimeSeries(key='your API key',output_format='pandas', indexing_type='date')

#read in ETF listing
ETF_df = dt.fread('ETF_symbols.csv')

#convert to pandas
ETF_df = ETF_df.to_pandas()

#loop through df
stock_price_history = pd.DataFrame()

for i in ETF_df.index:
    print(str(i) + " out of " + str(max(ETF_df.index)))
    try:
        stock = ETF_df.iloc[i,:][0]
        data, meta_data = ts.get_weekly_adjusted(symbol=stock)
        data.index.name = 'date'
        data.reset_index(inplace=True)
        data['symbol'] = stock
        stock_price_history = stock_price_history.append(data)
    except:
        pass
    
    
#write out csv
stock_price_history = dt.Frame(stock_price_history)

stock_price_history.to_csv('ETF_weekly_prices.csv')

n_year_length = 10 #this is how many years of returns we will analyze (we are interested in long term investments)
n_trading_weeks = 3 * 52 #we will only include ETFs that have traded for at least this many weeks (ideally, at least 3 years, or 156 weeks)

#read in ETF list that has already completed parsing
stocklist = fread('/home/kakasa09/stock_investing/ETF_weekly_prices.csv')

stocklist = stocklist %>% as_tibble() %>% select(date, `5. adjusted close`, `7. dividend amount`, symbol) %>%
  rename(Date = date,
         Adjusted = `5. adjusted close`,
         Dividend = `7. dividend amount`)

stocklist$Date <- as.Date(stocklist$Date)
latest_trading_date_global = max(stocklist$Date)

stocklist = stocklist %>% group_by(symbol) %>% mutate(trading_weeks = n_distinct(Date),
                                                      latest_trading_date = max(Date)) %>% ungroup() %>%
  filter(trading_weeks >= n_trading_weeks)  %>% filter(Date >= today() -
                                                         (365 * n_year_length)) %>%
  filter(latest_trading_date == latest_trading_date_global) %>% select(-latest_trading_date,-trading_weeks)

#this is how our raw data looks like after doing some filtering

kable(head(stocklist, n=5), format ='markdown')

#get the weekly log returns of our adjusted ETF prices    
#function modified from original source: https://www.business-science.io/investments/2016/11/30/Russell2000_Analysis.html

get_log_returns <- function(x, return_format = "tibble", period = 'weekly', ...) {
          # Convert tibble to xts
          x <- as.xts(x[, names(x) != "Date"], order.by = x$Date)
          # Get ETF prices
          log_returns_xts <- periodReturn(x = x$Adjusted, type = 'log', period = period, ...)
          # Rename
          names(log_returns_xts) <- "Log_Returns"
          stock_dates = index(log_returns_xts)
          # Return in xts format if tibble is not specified
     
            log_returns <- log_returns_xts %>%
              tibble::as_tibble() %>%
              mutate(Date = stock_dates) %>%
              mutate(Date = ymd(Date))

          log_returns
        }


# Use purrr to map the functions on the scraped ETF prices (scraping done in Python)
stocklist <- stocklist %>% group_by(symbol) %>% nest(.key = 'stock_prices') %>%
      mutate(
        log_returns  = map(stock_prices,
                           function(.x) get_log_returns(.x, return_format = "tibble")),
        mean_log_returns = map_dbl(log_returns, ~ mean(.$Log_Returns)),
        sd_log_returns   = map_dbl(log_returns, ~ sd(.$Log_Returns)),
        n_trade_weeks = map_dbl(stock_prices, nrow)
      )
    
# Calculate the average annual dividend payout:
dividend_profiles = stocklist %>% select(symbol, stock_prices) %>% unnest() %>% mutate(year  = year(Date)) %>% 
            group_by(symbol, year) %>% summarise(dividend_sum = sum(Dividend, na.rm = T)) %>%
              group_by(symbol) %>% summarise(dividend_payout_by_year_mean = mean(dividend_sum, na.rm = T)) %>% ungroup()
   

dividend_profiles[is.na(dividend_profiles)] <- 0

# Let's see how the first few rows look like:

kable(head(dividend_profiles, n=5), format = 'markdown')

#An ETF with high average log returns and low standard deviation of log returns will get a very high reward metric score.
#This makes intuitive sense: we are interested in ETFs with high returns that are also stable over time
#conversely, an ETF with a low average log return and high standard deviation will get a very low reward metric score

stocklist <- stocklist %>%
  mutate(reward_metric = (mean_log_returns * n_trade_weeks) / sd_log_returns) %>%
          left_join(dividend_profiles)


#Visualize the results with Plotly
#the darker the green, the better the reward metric is for the ETF

# Inputs
trade_day_thresh <- min(stocklist$n_trade_weeks)
lab <- "Select ETFs"
back_col <- '#2C3E50'
font_col <- '#FFFFFF'
line_col <- "#FFFFFF"
grid_col <- 'rgb(255, 255, 255)'
col_brew_pal <- 'YlGn'
# Plotly
plot_ly(data   = stocklist %>% filter(n_trade_weeks >= trade_day_thresh),
        type   = "scatter",
        mode   = "markers",
        x      = ~ sd_log_returns,
        y      = ~ mean_log_returns,
        color  = ~ reward_metric,
        colors = col_brew_pal,
        size   = ~ reward_metric,
        text   = ~ str_c("Ticker: ", symbol, "<br>",
                         "No. of Trading Weeks: ", n_trade_weeks, "<br>",
                         "Reward to Risk: ", round(reward_metric, 1)),
        marker = list(opacity = 0.8,
                      symbol = 'circle',
                      sizemode = 'diameter',
                      sizeref = 4.0,
                      line = list(width = 2, color = line_col))
) %>%
  layout(title   = str_c(lab, 'Analysis: ETF Risk vs Reward', sep = " "),
         xaxis   = list(title = 'Risk: StDev of Weekly Log Returns (SDWLR)',
                        gridcolor = grid_col,
                        zerolinewidth = 1,
                        ticklen = 5,
                        gridwidth = 2),
         yaxis   = list(title = 'Reward: Mean Weekly Log Returns (MWLR)',
                        gridcolor = grid_col,
                        zerolinewidth = 1,
                        ticklen = 5,
                        gridwith = 2),
         margin = list(l = 100,
                       t = 100,
                       b = 100),
         font   = list(color = font_col),
         paper_bgcolor = back_col,
         plot_bgcolor = back_col)

#limit time series to evaluating performance for the last 156 weeks (you can change this)
all_stock_cordata <- stocklist %>% 
  select(symbol, stock_prices) %>%
  unnest() %>% filter(Date >= today()-(n_trading_weeks*7))

#eliminate dates that are sparse
num_stocklist = length(unique(all_stock_cordata$symbol))

dates_sparsity = all_stock_cordata %>% group_by(Date) %>%
                          summarise(sparsity = 1- n_distinct(symbol)/num_stocklist ) %>% ungroup() %>%
                          filter(sparsity <= 0.01)

stocks_sparsity  = all_stock_cordata %>% filter(Date %in% unique(dates_sparsity$Date)) %>%
                                  group_by(symbol) %>% summarise(num_weeks = n()) %>%
                                    ungroup() %>% filter(num_weeks == dim(dates_sparsity)[1])

all_stock_cordata = all_stock_cordata %>% filter(symbol %in% unique(stocks_sparsity$symbol))

all_stock_cordata$adjusted_normalized = as.numeric(scale(  all_stock_cordata$Adjusted, scale =TRUE, center = TRUE))


# Spread format conducive to cor()
all_stock_cordata <- all_stock_cordata %>%
  select(symbol, Date, adjusted_normalized) %>% 
  spread(key = symbol, value = adjusted_normalized) %>%
  na.omit()

all_stock_cordata_matrix = t(as.matrix(all_stock_cordata[,-1]))

#Create 15 time series clusters: 15 seems to provide nice separation (be sure to try a few clusters to visually confirm)
#Clusters group ETFs together that have similar price trends over the last 3 years (they can also be a proxy for industry concentration within the ETFs)
num_ts_clust = 15

dtw_cluster = tsclust(all_stock_cordata_matrix, type="partitional",k=num_ts_clust,
                      distance="dtw_basic",centroid = "pam",seed=1234,trace=T,
                      args = tsclust_args(dist = list(window.size = 20)))

plot(dtw_cluster)

#visualize the slope of centroids with ggplot
clust_centroids = tibble(centroids = dtw_cluster@centroids) %>% mutate(cluster = 1:num_ts_clust) %>%
            unnest() %>% mutate(week_index = rep(1:(n()/num_ts_clust),num_ts_clust))


clust_centroids %>%
  ggplot(aes(x = week_index, y = centroids)) +
  geom_ref_line(h = 0) +
  geom_line(aes(col = cluster)) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~cluster, scales = "free_y") +
  theme(legend.position = "None", axis.text.x = element_text(angle=90)) +
  labs(title = "Cluster centroid trends over the last 3 years",
       x = "Week Index",
       y = "Normalized ETF Price")

# Function to return linear model of cluster centroids
centroid_model <- function(data) {
  lm(centroids ~ week_index, data = data)
}


# Function to return slope of the linear model
slope_centroid <- function(centroid_model) {
  centroid_model$coefficients[[2]]
}

cluster_slopes = clust_centroids %>% 
                group_by(cluster) %>% nest() %>%
                    mutate(
                      centroid_model = map(data, centroid_model),
                      slope = map_dbl(centroid_model, slope_centroid)) 

median_cluster_slope = median(cluster_slopes$slope, na.rm = T)

clusters_keep =   cluster_slopes %>% filter(slope > median_cluster_slope) %>%
                                                filter(slope > 0)

clusters_keep = unique(clusters_keep$cluster)

clusters_exclude = unique(clust_centroids$cluster)[unique(clust_centroids$cluster) %ni% clusters_keep]


#filter stocklist to only include the ETFs that belong to the clusters with positive trends
cluster_info = data.frame(symbol = names(dtw_cluster@datalist), ts_cluster = dtw_cluster@cluster) %>% 
                              filter(ts_cluster %in% clusters_keep)
  
stocklist = stocklist %>% inner_join(cluster_info)

top_n_limit <- 10  #select this many top performing ETFs within each time series cluster based on their Reward Metric ranking

hp <- stocklist %>% group_by(ts_cluster) %>% 
  mutate(rank = reward_metric %>% desc() %>% min_rank()) %>% ungroup() %>% 
  filter(rank <= top_n_limit) %>%
  arrange(rank)


# Function to return MWLR by year
means_by_year <- function(log_returns) {
  log_returns %>%
    mutate(year = year(Date)) %>%
    group_by(year) %>%
    summarize(mean_log_returns = mean(Log_Returns))
}


# Map function to data frame
hp <- hp %>%
  mutate(means_by_year = map(log_returns, means_by_year))


# Compute three attributes of high performing stocks and further filter---------------------------

# Attribute 1: Number of Times MWLR by Year Drops Below Zero
# Function to return number of times a stock's MDLR by year drops below zero
means_below_zero <- function(means_by_year) {
  means_by_year %>%
    filter(mean_log_returns < 0) %>%
    nrow()
}

# Map function to data frame for all stocks
hp <- hp %>%
  mutate(means_below_zero = map_dbl(means_by_year, means_below_zero))


# Attribute 2: Slope of MWLR by Year
# Function to return linear model
means_by_year_model <- function(means_by_year) {
  lm(mean_log_returns ~ year, data = means_by_year)
}


# Function to return slope of linear model
slope <- function(means_by_year_model) {
  means_by_year_model$coefficients[[2]]
}  


# Map modeling and slope functions
hp <- hp %>%
  mutate(
    means_by_year_model = map(means_by_year, means_by_year_model),
    slope = map_dbl(means_by_year_model, slope)
  )

# Attribute 3: Standard deviation of MWLR by year
sd_of_means_by_year <- function(means_by_year) {
  sd(means_by_year$mean_log_returns)
}

# Map to data frame
hp <- hp %>%
  mutate(sd_of_means_by_year = map_dbl(means_by_year, sd_of_means_by_year))


# Develop Growth-to-Consistency Metric -----------------------------------------

hp <- hp %>%
  mutate(growth_metric = slope /((means_below_zero + 1) * sd_of_means_by_year)) 


# Visualize Performance of Top N number of ETFs by Cluster --------------------------------------

top_n_limit <- 5

#Further trim list based on "Growth to Consistency" metric: select top N number of ETFs within each cluster

top_n_stock_profiles = hp %>% group_by(ts_cluster) %>% 
  mutate(rank = growth_metric %>% desc() %>% min_rank()) %>% ungroup() %>% 
  filter(rank <= top_n_limit) %>% 
  arrange(rank)

hp_stock_list = unique(top_n_stock_profiles$symbol)


ncol_hp_gg = 4
nrow_hp_gg = 5
num_pages_ggplot = ceiling(dim(top_n_stock_profiles)[1]/(ncol_hp_gg*nrow_hp_gg))

for(i in 1:num_pages_ggplot){
  x= top_n_stock_profiles %>%
    mutate(symbol_cluster = paste("(",ts_cluster,") ",symbol, sep = "")) %>% 
    select(ts_cluster, symbol_cluster, means_by_year)  %>% 
    unnest() %>%
    mutate(ts_cluster = as.factor(ts_cluster)) %>% rename(ETF_cluster = ts_cluster) %>% 
    ggplot(aes(x = year, y = mean_log_returns)) +
    geom_ref_line(h = 0) +
    geom_line(aes(col = ETF_cluster)) +
    geom_smooth(method = "lm", se = FALSE) +
    facet_wrap_paginate(facets = 'symbol_cluster', ncol = ncol_hp_gg, nrow = nrow_hp_gg, page = i, scales = "free_y") +
    theme(legend.position = "None", axis.text.x = element_text(angle=90)) +
    labs(title = paste0("Average MWLR per annum for select list of ",top_n_limit*length(clusters_keep), " (top ",top_n_limit," ETFs per cluster)"),
         subtitle = paste0("page ",i," of ",num_pages_ggplot),
         x = "Year",
         y = "Avg. Annual Mean Weekly Log Returns (MWLR)")
  print(x) }

top_n_stock_profiles %>%
  mutate(symbol_cluster = paste("(",ts_cluster,") ",symbol, sep = "")) %>% 
  select(ts_cluster, symbol_cluster, stock_prices)  %>% 
  unnest() %>%
  mutate(ts_cluster = as.factor(ts_cluster)) %>% rename(ETF_cluster = ts_cluster) %>% 
  ggplot(aes(x = Date, y = Adjusted, col = ETF_cluster)) +
  geom_line() +
  facet_wrap(~ symbol_cluster, nrow = 5, scales = "free_y") +
  labs(title = paste0("Prospecting Best ETFs from selected list of ",top_n_limit*length(clusters_keep), " (top ",top_n_limit," per cluster)"),
       subtitle = paste0(top_n_limit, " ETFs with best Growth, Most Consistency within each of our focus clusters"),
       x = "Year",
       y = "Price per Share")

#apply weights (we will weight dividend performance at 5%)
dividend_wt = 0.05
stock_growth_wt = 1-dividend_wt

#re-rank it, but without the clusters this time
top_n_stock_profiles2 = hp %>% filter(symbol %in% hp_stock_list) %>%
  mutate(rank = growth_metric %>% desc() %>% min_rank()) %>% 
  select(ts_cluster, symbol, dividend_payout_by_year_mean, rank) %>%
  arrange(rank) %>% rename(stock_growth_rank = rank) %>%
                        mutate(dividend_payout_rank = row_number(desc(dividend_payout_by_year_mean))) %>%
                        mutate(final_rank = stock_growth_wt*stock_growth_rank+
                                 dividend_payout_rank * dividend_wt) %>%
                          mutate(final_rank = row_number(final_rank)) %>%
                          arrange(final_rank) %>% select(ts_cluster, symbol, final_rank, stock_growth_rank, dividend_payout_rank)



# Purchase Qty Recommendation --------------------------------------
top_n_recent_prices = hp %>% filter(symbol %in% hp_stock_list) %>%
                            select(symbol, stock_prices) %>% unnest() 

max_date= max(top_n_recent_prices$Date)

top_n_recent_prices = top_n_recent_prices %>% filter(Date == max_date) %>%
                      select(symbol, Adjusted) %>% rename(recent_price = Adjusted)


#specify how many ETFs to consider 
n_stocks_consider = 10

#Specify # of ETFs we want to invest in:
final_invest_ETF_count = 5

#recommendation on what ETFs to consider: 
top_n_stock_profiles_consider = top_n_stock_profiles2 %>% left_join(top_n_recent_prices) %>% filter(final_rank <= n_stocks_consider)

consideration_stock_list = unique(top_n_stock_profiles_consider$symbol)

kable(top_n_stock_profiles_consider, format = 'markdown')

top_n_stock_profiles %>% filter(symbol %in% consideration_stock_list) %>%
  select(symbol, log_returns) %>%
  unnest() %>% spread(key = symbol, value = Log_Returns) %>%
  na.omit()  %>% select(-Date) %>%
  cor() %>%
  corrplot(order   = "hclust", method = 'number',
           addrect = final_invest_ETF_count, rect.col = 'yellow', rect.lwd = 5, number.digits = 2)

#final list based on our analysis
final_invest_stock_list = c('IHI','VGT','IGV','XLG','VPU')


#make any manual additions for specific ETFs you are interested in (in this case, I'm interested in water ETFs) 
ETF_additions = data.frame(symbol = c('CGW', 'TBLU'),
                           recent_price = c( 37.46, 31), stringsAsFactors = F)


final_invest_stock_list = c(final_invest_stock_list, ETF_additions$symbol)



#Correlation plot of final ETF selections:
final_cor_data = fread('/home/kakasa09/stock_investing/ETF_weekly_prices.csv')

final_cor_data = final_cor_data %>% as_tibble() %>% select(date, `5. adjusted close`, `7. dividend amount`, symbol) %>%
  rename(Date = date,
         Adjusted = `5. adjusted close`,
         Dividend = `7. dividend amount`)

final_cor_data$Date <- as.Date(final_cor_data$Date )    

final_cor_data = final_cor_data %>%  filter(Date >= today()-(365*n_year_length)) %>%
                  filter(symbol %in% c(final_invest_stock_list)) %>% group_by(symbol) %>% nest(.key = 'stock_prices') %>%
                    mutate(log_returns  = map(stock_prices,
                                              function(.x)
                                                get_log_returns(.x, return_format = "tibble"))) 

final_cor_data %>% select(symbol, log_returns) %>%
                        unnest() %>% spread(key = symbol, value = Log_Returns) %>%
                          na.omit()  %>% select(-Date) %>%
                              cor() %>%
                                  corrplot(order = "hclust", method = 'number', number.digits = 2)

#visualize price trends for our final ETF list
final_cor_data %>% 
  mutate(symbol = if_else(symbol %in% ETF_additions$symbol, paste0(symbol, " (manual)"),
                          paste0(symbol, " (auto)"))) %>% 
  select(symbol, stock_prices)  %>% 
  unnest() %>%
  ggplot(aes(x = Date, y = Adjusted, col = symbol)) +
  geom_line() +
  facet_wrap(~ symbol,  scales = "free") +
  labs(title = paste0("Price trends for final list of ETFs"),
       subtitle = paste0("Weekly price/share trends for the past ", n_year_length, " years"),
       x = "Year",
       y = "Price per Share") + theme(legend.position = 'none')