Install and call required packages

library(readxl)
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(corrplot)
## corrplot 0.92 loaded
library(elasticnet)
## Loading required package: lars
## Loaded lars 1.3

Read stocks and fetch data

table <- read_excel("stock_data_outputV3.xlsx")
## New names:
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
stock_list <- table[[1]]
stock_len <- length(stock_list)
df_price <- getSymbols(stock_list[1],auto.assign = F,from = "2023-10-02",to = "2024-03-28")[,4]
for (i in 2:stock_len){
  df_price <- cbind(df_price,getSymbols(stock_list[i],auto.assign = F,from = "2023-10-02",to = "2024-03-28")[,4])
}
colnames(df_price) <- stock_list
df_price <- data.frame(df_price)
head(df_price)
##              NVDA PAYS   APP   ASML   AVGO  ACMR  STNE    ODD   ANET    CRM
## 2023-10-02 447.82 1.86 40.82 582.18 834.95 18.62 10.68 28.080 187.81 203.71
## 2023-10-03 435.17 1.83 39.13 570.60 814.83 17.15 10.25 29.570 184.89 199.83
## 2023-10-04 440.41 1.83 40.33 590.86 823.81 18.00 10.50 27.100 186.03 201.87
## 2023-10-05 446.88 1.83 39.70 581.69 823.86 17.30 10.08 25.920 188.03 202.01
## 2023-10-06 457.62 1.85 40.38 597.36 845.26 18.61 10.01 29.645 194.02 207.36
## 2023-10-09 452.73 1.84 40.85 591.37 855.88 18.22  9.91 29.700 196.40 207.22
##              FICO   GRMN    NOW    PTC   CDNS   WDAY  AGYS    DV   MANH  DSGX
## 2023-10-02 863.59 103.51 555.23 140.84 236.49 213.66 66.13 27.88 198.86 72.74
## 2023-10-03 836.54 103.51 541.24 137.84 228.69 207.30 66.02 27.38 195.13 71.93
## 2023-10-04 860.85 104.10 548.64 139.63 235.05 209.37 67.13 27.89 199.34 71.75
## 2023-10-05 862.53 104.59 547.17 139.82 234.77 205.14 66.99 27.52 199.51 72.21
## 2023-10-06 873.52 106.13 560.72 140.34 243.95 213.63 70.85 28.03 206.85 74.06
## 2023-10-09 889.39 107.17 564.43 139.94 247.34 217.75 71.48 28.13 208.21 74.68
##              INTU   ADSK   DJCO   SNPS   KARO   PANW   ANSS  FTNT  EXLS    BMI
## 2023-10-02 517.80 208.99 293.39 463.82 20.395 236.77 297.40 58.63 28.03 142.91
## 2023-10-03 507.48 203.28 291.21 450.72 18.590 231.71 287.39 57.61 27.51 142.30
## 2023-10-04 517.46 207.06 292.00 461.32 19.940 236.44 290.82 58.26 27.67 142.48
## 2023-10-05 516.30 203.96 292.50 459.75 21.050 236.78 290.18 58.02 27.67 142.36
## 2023-10-06 533.30 210.04 296.29 473.26 21.580 246.69 296.51 59.93 28.44 142.99
## 2023-10-09 534.91 210.56 297.00 479.92 21.690 252.78 300.91 59.51 29.22 143.03

Calculate returns and see correlation plot

df_daily_return <- (df_price[2:dim(df_price)[1],] - df_price[1:dim(df_price)[1]-1,])/df_price[1:dim(df_price)[1]-1,]
corrplot(cor(df_daily_return),method = "color")

PCA and analysis

pca_stock <- prcomp(df_daily_return,scale. = T)
biplot(pca_stock)

var_explained <- pca_stock$sdev**2/sum(pca_stock$sdev**2)
plot(1:30,var_explained,type = "o")
text(c(3,4),var_explained[1:2],labels = round(var_explained[1:2],3))

cumsum(pca_stock$sdev**2/sum(pca_stock$sdev**2))
##  [1] 0.3422001 0.4078920 0.4567652 0.5014547 0.5433820 0.5828582 0.6189820
##  [8] 0.6540478 0.6856207 0.7161329 0.7451372 0.7713715 0.7962241 0.8184959
## [15] 0.8390167 0.8578738 0.8761883 0.8920298 0.9065093 0.9204045 0.9331354
## [22] 0.9442869 0.9545145 0.9643615 0.9728734 0.9810923 0.9868222 0.9921839
## [29] 0.9971099 1.0000000
which(cumsum(pca_stock$sdev**2/sum(pca_stock$sdev**2))>0.9)
##  [1] 19 20 21 22 23 24 25 26 27 28 29 30

Export first 19 PCs explaining 90% variance (can choose fewer PCs)

pca_coords <- data.frame(pca_stock$x[,1:19])
write.csv(pca_coords,file = "pca_coords.csv")

Sparse PCA

  1. The idea here is to filter out variables that are important when doing PCA analysis;
  2. When doing ordinary PCA, we need 19 principal components to cover 90% of total variance so 19 components is a good start to find SPCA as variance explained increases with increasing loadings and components;
n_component <- seq(19,30)
n_loading <- seq(1,30)
spca_test_var_explained_vec <- c()
stop = FALSE
for(num_coponent in n_component){
  for (num_loading in n_loading){
    spca_stock_test <- spca(df_daily_return,num_coponent,para = rep(num_loading,num_coponent),
                            type = "predictor",sparse = "varnum")
    spca_stock_test_var_explained <- sum(spca_stock_test$pev)
    spca_test_var_explained_vec <- append(spca_test_var_explained_vec,spca_stock_test_var_explained)
    length_vec <- length(spca_test_var_explained_vec)
    if (spca_stock_test_var_explained >0.9){
      stop = TRUE
      break
    }
  }
  if (stop){break}
}
print(length(spca_test_var_explained_vec))#length is 25. So n_component = 19, n_loading =25
## [1] 25
spca_stock <- spca(df_daily_return,19,para = rep(25,19),type = "predictor",sparse = "varnum")
print(sum(spca_stock$pev))
## [1] 0.9030985
  1. Get the principal components with some companies’ loadings to be zero
spca_coods <- as.matrix(df_daily_return)%*%spca_stock$loadings
write.csv(spca_coods,"spca_coods.csv")