library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.1.0 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(ggplot2)
library(csv)
#1. Import the data
stock <- read.csv('/Users/macbookair/Downloads/m-fac9003.csv')
head(stock)
## AA AGE CAT F FDX GM HPQ KMB MEL NYT PG
## 1 -16.40 -12.17 -4.44 -0.06 -2.28 -2.12 -6.19 -11.01 -10.77 -6.30 -8.89
## 2 4.04 4.95 8.84 6.02 10.47 8.97 -4.01 -5.20 0.34 -4.62 -0.84
## 3 0.12 13.08 0.17 2.06 10.84 1.57 5.67 3.21 -0.17 -0.66 5.41
## 4 -4.28 -11.06 0.25 -5.67 -2.44 -4.19 -5.29 -0.65 -2.20 -10.60 4.26
## 5 5.81 19.70 8.52 3.89 -16.17 10.94 8.81 8.83 11.85 11.59 16.35
## 6 -4.05 -1.44 -22.10 -5.79 -2.81 -2.70 -1.47 1.55 -7.76 -0.12 4.80
## TRB TXN SP500
## 1 -13.04 -7.61 -7.52
## 2 -0.37 4.97 0.21
## 3 2.36 2.69 1.77
## 4 -7.98 -6.85 -3.34
## 5 8.82 22.88 8.55
## 6 -0.64 -5.87 -1.53
#2. Calculate the number of periods that is collected
num = dim(stock)[1]
num
## [1] 168
#3. Choose and Remove S&P500 market index monthly return from the data
market = stock[,14]
stock13 = stock[,c(-14)]
head(stock13)
## AA AGE CAT F FDX GM HPQ KMB MEL NYT PG
## 1 -16.40 -12.17 -4.44 -0.06 -2.28 -2.12 -6.19 -11.01 -10.77 -6.30 -8.89
## 2 4.04 4.95 8.84 6.02 10.47 8.97 -4.01 -5.20 0.34 -4.62 -0.84
## 3 0.12 13.08 0.17 2.06 10.84 1.57 5.67 3.21 -0.17 -0.66 5.41
## 4 -4.28 -11.06 0.25 -5.67 -2.44 -4.19 -5.29 -0.65 -2.20 -10.60 4.26
## 5 5.81 19.70 8.52 3.89 -16.17 10.94 8.81 8.83 11.85 11.59 16.35
## 6 -4.05 -1.44 -22.10 -5.79 -2.81 -2.70 -1.47 1.55 -7.76 -0.12 4.80
## TRB TXN
## 1 -13.04 -7.61
## 2 -0.37 4.97
## 3 2.36 2.69
## 4 -7.98 -6.85
## 5 8.82 22.88
## 6 -0.64 -5.87
#4. Convert data into martrix
stock13 = as.matrix(stock13)
#5: Find beta, e, D and R-square
n = dim(stock13)[2]
n
## [1] 13
cal = rep(1,num)
head(cal)
## [1] 1 1 1 1 1 1
X = cbind(cal, market)
head(X)
## cal market
## [1,] 1 -7.52
## [2,] 1 0.21
## [3,] 1 1.77
## [4,] 1 -3.34
## [5,] 1 8.55
## [6,] 1 -1.53
betta_hat = solve(t(X)%*%X)%*%t(X)%*%stock13
betta_hat
## AA AGE CAT F FDX GM HPQ
## cal 0.549124 0.7218061 0.8393521 0.4543643 0.7995790 0.1982025 0.6835681
## market 1.291591 1.5141359 0.9406928 1.2192453 0.8051166 1.0457019 1.6279512
## KMB MEL NYT PG TRB TXN
## cal 0.5463020 0.8849263 0.4904120 0.8880914 0.6512465 1.438887
## market 0.5498052 1.1228708 0.7706495 0.4688034 0.7178808 1.796412
E_hat = stock13 - X%*%betta_hat
head(E_hat)
## AA AGE CAT F FDX GM
## [1,] -7.2363588 -1.5055042 1.7946575 8.6543606 2.9748981 5.5454755
## [2,] 3.2196419 3.9102254 7.8031024 5.3095942 9.5013465 8.5522001
## [3,] -2.7152402 9.6781734 -2.3343784 -0.5524285 8.6153645 -0.4790948
## [4,] -0.5152096 -6.7245922 2.5525617 -2.0520849 -0.5504894 -0.8955583
## [5,] -5.7822280 6.0323321 -0.3622754 -6.9889119 -23.8533263 1.8010466
## [6,] -2.6229896 0.1548218 -21.5000922 -4.3789190 -2.3777506 -1.2982786
## HPQ KMB MEL NYT PG TRB TXN
## [1,] 5.3686247 -7.4217666 -3.2109382 -0.9951281 -6.252690 -8.2927829 4.460129
## [2,] -5.0354379 -5.8617611 -0.7807291 -5.2722484 -1.826540 -1.1720014 3.153867
## [3,] 2.1049583 1.6905428 -3.0424075 -2.5144615 3.692127 0.4381045 -1.928535
## [4,] -0.5362112 0.6400475 0.6654621 -8.5164428 4.937712 -6.2335246 -2.288872
## [5,] -5.7925506 3.5828633 1.3645288 4.5105352 11.453640 2.0308727 6.081793
## [6,] 0.3371972 1.8449000 -6.9269340 0.5686817 4.629178 -0.1928888 -4.560377
diagD_hat = diag(t(E_hat)%*%E_hat)/(num-2)
stockvar = apply(stock13, 2, var)
r_square = 1- diag(t(E_hat)%*%E_hat)/((num-1)*stockvar)
res_std = sqrt(diagD_hat)
cov_factor = var(market)*t(betta_hat)%*%betta_hat + diag(diagD_hat) #omega
#compute global minimum varience portfolio weight
one.vec = rep(1,13)
a = solve(cov_factor)%*%one.vec
b = t(one.vec)%*%a
mvpw = a / as.numeric(b)
mvpw
## [,1]
## AA 0.03501615
## AGE -0.03373670
## CAT 0.05314427
## F 0.05576422
## FDX 0.06451069
## GM 0.12369835
## HPQ -0.03290554
## KMB 0.28022778
## MEL 0.01901912
## NYT 0.19373120
## PG 0.19019732
## TRB 0.14314200
## TXN -0.09180887
#Plot the weight for each MVP
barplot(t(mvpw), horiz=F, main="Weight for each MVP", col="yellow", cex.names = 0.70, las=2)
