1. Download ETF daily data

  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

## Version 0.4-0 included new data defaults. See ?getSymbols.

  library(plyr)

  library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following objects are masked from 'package:xts':
## 
##     first, last

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

  library(PerformanceAnalytics)

## 
## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

  library(tidyverse)

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

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

## -- Conflicts ----------------------------------------------- tidyverse_conflicts() --
## x dplyr::arrange()   masks plyr::arrange()
## x purrr::compact()   masks plyr::compact()
## x dplyr::count()     masks plyr::count()
## x dplyr::failwith()  masks plyr::failwith()
## x dplyr::filter()    masks stats::filter()
## x dplyr::first()     masks xts::first()
## x dplyr::id()        masks plyr::id()
## x dplyr::lag()       masks stats::lag()
## x dplyr::last()      masks xts::last()
## x dplyr::mutate()    masks plyr::mutate()
## x dplyr::rename()    masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()

  library(tidyr)

  tickers <- c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD")

  n <- length(tickers)

  data <- new.env()

  getSymbols(tickers, src="yahoo", from = "2010-01-01",
           env = data, auto.assign = T)

## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.

## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols

## [1] "SPY" "QQQ" "EEM" "IWM" "EFA" "TLT" "IYR" "GLD"

  con = gzcon(url('http://www.systematicportfolio.com/sit.gz','rb'))

  source(con)

  close(con)

  bt.prep(data, align = 'remove.na', fill.gaps = T)

  prices <- data$prices

  ret <- prices/mlag(prices) - 1

  ret <- ret[-1, ]

2. Daily returns, covariance matrix and MVP

  ret.sample <- ret['2010-01-01/2010-12-30', ]

  avg.ret <- colMeans(coredata(ret.sample))

  cov.d <- cov(coredata(ret.sample))

  ones.vec <- matrix(rep(1,n), ncol = 1)

  numerator <- solve(cov.d) %*% ones.vec

  denominator <- t(ones.vec)%*%numerator

  mvp.d <- numerator / as.numeric(denominator)

mvp.d

##            [,1]
## EEM -0.10766006
## EFA -0.15032757
## GLD  0.10316006
## IWM -0.25701223
## IYR -0.08291293
## QQQ  0.03110772
## SPY  1.06528877
## TLT  0.39835624

3. Weekly returns, covariance matrix and MVP

  ret.sample.w <-   ret.sample %>% to.weekly(indexAt = "lastof", OHLC = FALSE)

  avg.ret.w <- colMeans(coredata(ret.sample.w))

  cov.d.w <- cov(coredata(ret.sample.w))

  numerator.w <- solve(cov.d.w) %*% ones.vec

  denominator.w <- t(ones.vec)%*%numerator.w

  mvp.w <- numerator.w / as.numeric(denominator.w)

  mvp.w

##            [,1]
## EEM -0.18728208
## EFA -0.03000348
## GLD  0.07266435
## IWM -0.26420190
## IYR -0.19689304
## QQQ -0.25396710
## SPY  1.46621438
## TLT  0.39346887

4. Monthly returns, covariance matrix and MVP

  ret.sample.m <-   ret.sample %>% to.monthly(indexAt = "lastof", OHLC = FALSE)

  avg.ret.m <- colMeans(coredata(ret.sample.m))

  cov.d.m <- cov(coredata(ret.sample.m))

  numerator.m <- solve(cov.d.m) %*% ones.vec

  denominator.m <- t(ones.vec)%*%numerator.m

  mvp.m <- numerator.m / as.numeric(denominator.m)

  mvp.m

##           [,1]
## EEM -0.5289397
## EFA -0.4680952
## GLD -0.4345513
## IWM -1.0475222
## IYR -0.2347427
## QQQ -0.1207826
## SPY  3.3755502
## TLT  0.4590836