#========================================
# Error message when running code as follow:
# Error in t(X) %*% as.matrix(X) : 
  # requires numeric/complex matrix/vector arguments
#========================================
rm(list=ls())
retdata <- read.csv("C:/Users/HP/Downloads/m-fac9003.csv", header=FALSE)
t = dim(retdata)[1]
t
## [1] 169
market = retdata[,14]

#Get 3-Month Treasury Bill Secondary Market Rate from the FRED
riskfree = retdata[,15] 

retdata1 = retdata[,c(-14, -15)]
retdata1 = as.matrix(retdata1)
n = dim(retdata1)[2]
n
## [1] 13
ones = rep(1,t)
X = cbind(ones,market)
as.matrix(X)
##        ones market  
##   [1,] "1"  "SP500" 
##   [2,] "1"  "-7.52" 
##   [3,] "1"  "0.21"  
##   [4,] "1"  "1.77"  
##   [5,] "1"  "-3.34" 
##   [6,] "1"  "8.55"  
##   [7,] "1"  "-1.53" 
##   [8,] "1"  "-1.16" 
##   [9,] "1"  "-10.05"
##  [10,] "1"  "-5.73" 
##  [11,] "1"  "-1.27" 
##  [12,] "1"  "5.4"   
##  [13,] "1"  "1.92"  
##  [14,] "1"  "3.63"  
##  [15,] "1"  "6.23"  
##  [16,] "1"  "1.73"  
##  [17,] "1"  "-0.44" 
##  [18,] "1"  "3.4"   
##  [19,] "1"  "-5.25" 
##  [20,] "1"  "4.02"  
##  [21,] "1"  "1.52"  
##  [22,] "1"  "-2.35" 
##  [23,] "1"  "0.77"  
##  [24,] "1"  "-4.77" 
##  [25,] "1"  "10.82" 
##  [26,] "1"  "-2.31" 
##  [27,] "1"  "0.64"  
##  [28,] "1"  "-2.52" 
##  [29,] "1"  "2.48"  
##  [30,] "1"  "-0.21" 
##  [31,] "1"  "-2.04" 
##  [32,] "1"  "3.67"  
##  [33,] "1"  "-2.66" 
##  [34,] "1"  "0.67"  
##  [35,] "1"  "-0.03" 
##  [36,] "1"  "2.77"  
##  [37,] "1"  "0.74"  
##  [38,] "1"  "0.46"  
##  [39,] "1"  "0.8"   
##  [40,] "1"  "1.62"  
##  [41,] "1"  "-2.78" 
##  [42,] "1"  "2.03"  
##  [43,] "1"  "-0.18" 
##  [44,] "1"  "-0.79" 
##  [45,] "1"  "3.19"  
##  [46,] "1"  "-1.24" 
##  [47,] "1"  "1.69"  
##  [48,] "1"  "-1.55" 
##  [49,] "1"  "0.75"  
##  [50,] "1"  "3"     
##  [51,] "1"  "-3.28" 
##  [52,] "1"  "-4.87" 
##  [53,] "1"  "0.85"  
##  [54,] "1"  "0.9"   
##  [55,] "1"  "-3.03" 
##  [56,] "1"  "2.79"  
##  [57,] "1"  "3.39"  
##  [58,] "1"  "-3.08" 
##  [59,] "1"  "1.67"  
##  [60,] "1"  "-4.39" 
##  [61,] "1"  "0.76"  
##  [62,] "1"  "1.95"  
##  [63,] "1"  "3.13"  
##  [64,] "1"  "2.26"  
##  [65,] "1"  "2.33"  
##  [66,] "1"  "3.16"  
##  [67,] "1"  "1.67"  
##  [68,] "1"  "2.73"  
##  [69,] "1"  "-0.48" 
##  [70,] "1"  "3.57"  
##  [71,] "1"  "-0.94" 
##  [72,] "1"  "3.66"  
##  [73,] "1"  "1.32"  
##  [74,] "1"  "2.85"  
##  [75,] "1"  "0.29"  
##  [76,] "1"  "0.38"  
##  [77,] "1"  "0.93"  
##  [78,] "1"  "1.87"  
##  [79,] "1"  "-0.2"  
##  [80,] "1"  "-5"    
##  [81,] "1"  "1.46"  
##  [82,] "1"  "4.99"  
##  [83,] "1"  "2.2"   
##  [84,] "1"  "6.92"  
##  [85,] "1"  "-2.56" 
##  [86,] "1"  "5.71"  
##  [87,] "1"  "0.18"  
##  [88,] "1"  "-4.69" 
##  [89,] "1"  "5.41"  
##  [90,] "1"  "5.44"  
##  [91,] "1"  "3.93"  
##  [92,] "1"  "7.39"  
##  [93,] "1"  "-6.17" 
##  [94,] "1"  "4.9"   
##  [95,] "1"  "-3.86" 
##  [96,] "1"  "4.03"  
##  [97,] "1"  "1.14"  
##  [98,] "1"  "0.6"   
##  [99,] "1"  "6.62"  
## [100,] "1"  "4.58"  
## [101,] "1"  "0.5"   
## [102,] "1"  "-2.3"  
## [103,] "1"  "3.53"  
## [104,] "1"  "-1.58" 
## [105,] "1"  "-14.99"
## [106,] "1"  "5.86"  
## [107,] "1"  "7.7"   
## [108,] "1"  "5.55"  
## [109,] "1"  "5.27"  
## [110,] "1"  "3.74"  
## [111,] "1"  "-3.6"  
## [112,] "1"  "3.51"  
## [113,] "1"  "3.44"  
## [114,] "1"  "-2.87" 
## [115,] "1"  "5.06"  
## [116,] "1"  "-3.58" 
## [117,] "1"  "-1.02" 
## [118,] "1"  "-3.24" 
## [119,] "1"  "5.85"  
## [120,] "1"  "1.48"  
## [121,] "1"  "5.35"  
## [122,] "1"  "-5.53" 
## [123,] "1"  "-2.47" 
## [124,] "1"  "9.2"   
## [125,] "1"  "-3.55" 
## [126,] "1"  "-2.67" 
## [127,] "1"  "1.92"  
## [128,] "1"  "-2.13" 
## [129,] "1"  "5.56"  
## [130,] "1"  "-5.85" 
## [131,] "1"  "-1"    
## [132,] "1"  "-8.52" 
## [133,] "1"  "-0.08" 
## [134,] "1"  "3.03"  
## [135,] "1"  "-9.64" 
## [136,] "1"  "-6.79" 
## [137,] "1"  "7.36"  
## [138,] "1"  "0.21"  
## [139,] "1"  "-2.79" 
## [140,] "1"  "-1.37" 
## [141,] "1"  "-6.69" 
## [142,] "1"  "-8.39" 
## [143,] "1"  "1.63"  
## [144,] "1"  "7.36"  
## [145,] "1"  "0.62"  
## [146,] "1"  "-1.69" 
## [147,] "1"  "-2.22" 
## [148,] "1"  "3.52"  
## [149,] "1"  "-6.29" 
## [150,] "1"  "-1.05" 
## [151,] "1"  "-7.39" 
## [152,] "1"  "-8.04" 
## [153,] "1"  "0.35"  
## [154,] "1"  "-11.14"
## [155,] "1"  "8.51"  
## [156,] "1"  "5.6"   
## [157,] "1"  "-6.13" 
## [158,] "1"  "-2.84" 
## [159,] "1"  "-1.8"  
## [160,] "1"  "0.74"  
## [161,] "1"  "8.01"  
## [162,] "1"  "5"     
## [163,] "1"  "1.06"  
## [164,] "1"  "1.55"  
## [165,] "1"  "1.71"  
## [166,] "1"  "-1.27" 
## [167,] "1"  "5.42"  
## [168,] "1"  "0.64"  
## [169,] "1"  "5"
#b_hat = solve(t(X)%*% as.matrix(X)) %*% t(X) %*% retdata1
#E_hat = retdata1 - X%*%b_hat
#diagD_hat = diag(t(E_hat)%*%E_hat)/(t-2)

#========================================
#R-square ----
#========================================
#retvar = apply(retdata1,2,var) 
#R_2 = 1 - diag(t(E_hat)%*%E_hat)/((t-1)*retvar)
#res_std = sqrt(diagD_hat)
#cov_factor = var(market)*t(b_hat)%*%b_hat + diag(diagD_hat) 
#sd = sqrt(diag(cov_factor));
#cor_factor = cov_factor/(sd%*%t(sd));

#========================================
# sample variance and correlation matrix
#cov_sample = cov(retdata1);
#cor_sample = cor(retdata1);

#========================================
# use factor covariance matrix to compute global minimum variance portfolio
#========================================
#one.vec = rep(1,15)
#a = solve(cov_factor)%*%one.vec
#b = t(one.vec)%*%a
#mvp.w =a / as.numeric(b)
#as.vector(mvp.w, names = rownames(mvp.w))
#barplot(as.vector(mvp.w), ylim = c(-0.01, 0.4), names.arg = rownames(mvp.w), cex.names = 0.5)
# ggplot2
library(ggplot2)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble  3.1.5     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.0.2     v forcats 0.5.1
## v purrr   0.3.4
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
#
#tickers<-rownames(mvp.w)
#weighti<-as.numeric(mvp.w)
#mvp.w.df<-data.frame(tickers, weighti)
#
#mvp.w.df %>% 
  #ggplot(aes(tickers, weighti)) +
  #geom_bar(stat = "identity")