rm(list=ls())

library(fBasics)

## Loading required package: timeDate

## Loading required package: timeSeries

library(tidyverse)

## -- Attaching packages ------------------------- tidyverse 1.3.0 --

## v ggplot2 3.3.0     v purrr   0.3.3
## v tibble  2.1.3     v dplyr   0.8.3
## v tidyr   1.0.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0

## -- Conflicts ---------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks timeSeries::filter(), stats::filter()
## x dplyr::lag()    masks timeSeries::lag(), stats::lag()

library(csv)

retdata = read.csv('m-fac9003.csv')

t = dim(retdata)[1]

market = retdata[,14]

str(retdata)

## 'data.frame':    168 obs. of  14 variables:
##  $ AA   : num  -16.4 4.04 0.12 -4.28 5.81 ...
##  $ AGE  : num  -12.17 4.95 13.08 -11.06 19.7 ...
##  $ CAT  : num  -4.44 8.84 0.17 0.25 8.52 ...
##  $ F    : num  -0.06 6.02 2.06 -5.67 3.89 ...
##  $ FDX  : num  -2.28 10.47 10.84 -2.44 -16.17 ...
##  $ GM   : num  -2.12 8.97 1.57 -4.19 10.94 ...
##  $ HPQ  : num  -6.19 -4.01 5.67 -5.29 8.81 ...
##  $ KMB  : num  -11.01 -5.2 3.21 -0.65 8.83 ...
##  $ MEL  : num  -10.77 0.34 -0.17 -2.2 11.85 ...
##  $ NYT  : num  -6.3 -4.62 -0.66 -10.6 11.59 ...
##  $ PG   : num  -8.89 -0.84 5.41 4.26 16.35 ...
##  $ TRB  : num  -13.04 -0.37 2.36 -7.98 8.82 ...
##  $ TXN  : num  -7.61 4.97 2.69 -6.85 22.88 ...
##  $ SP500: num  -7.52 0.21 1.77 -3.34 8.55 ...

retdata$date

## NULL

riskfree = retdata[,12]

retdata1 = retdata[,c(-14, -12)]

retdata1 = as.matrix(retdata1)

n = dim(retdata1)[2]

ones = rep(1,t)

X = cbind(ones,market)

b_hat = solve(t(X)%*%X)%*%t(X)%*%retdata1

E_hat = retdata1 - X%*%b_hat

diagD_hat = diag(t(E_hat)%*%E_hat)/(t-2)

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));

cov_sample = cov(retdata1);

cor_sample = cor(retdata1);

one.vec = rep(1,12)

a = solve(cov_factor)%*%one.vec

b = t(one.vec)%*%a

mvp.w = a/as.numeric(b)

mvp.w

##            [,1]
## AA   0.04412460
## AGE -0.02676160
## CAT  0.06481505
## F    0.06427293
## FDX  0.07424653
## GM   0.13385950
## HPQ -0.02875738
## KMB  0.30455158
## MEL  0.03522196
## NYT  0.21197141
## PG   0.21228558
## TXN -0.08983016