con = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)
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
## 
## Attaching package: 'TTR'
## The following object is masked _by_ '.GlobalEnv':
## 
##     DVI
## Version 0.4-0 included new data defaults. See ?getSymbols.
#
firm_data1 = read.csv('3firmExample_data3.csv')
str(firm_data1)
## 'data.frame':    59 obs. of  4 variables:
##  $ date     : Factor w/ 59 levels "1995/10/1","1995/11/1",..: 4 5 6 7 8 9 10 1 2 3 ...
##  $ Nordstrom: num  -0.03615 -0.0568 0.07821 -0.00302 -0.02757 ...
##  $ Starbucks: num  0.00521 -0.02105 0.21244 0.2036 0.04797 ...
##  $ Microsoft: num  0.1213 0.13923 0.03529 0.06501 0.00138 ...
firm_data1$date
##  [1] 1995/3/1  1995/4/1  1995/5/1  1995/6/1  1995/7/1  1995/8/1  1995/9/1 
##  [8] 1995/10/1 1995/11/1 1995/12/1 1996/1/1  1996/2/1  1996/3/1  1996/4/1 
## [15] 1996/5/1  1996/6/1  1996/7/1  1996/8/1  1996/9/1  1996/10/1 1996/11/1
## [22] 1996/12/1 1997/1/1  1997/2/1  1997/3/1  1997/4/1  1997/5/1  1997/6/1 
## [29] 1997/7/1  1997/8/1  1997/9/1  1997/10/1 1997/11/1 1997/12/1 1998/1/1 
## [36] 1998/2/1  1998/3/1  1998/4/1  1998/5/1  1998/6/1  1998/7/1  1998/8/1 
## [43] 1998/9/1  1998/10/1 1998/11/1 1998/12/1 1999/1/1  1999/2/1  1999/3/1 
## [50] 1999/4/1  1999/5/1  1999/6/1  1999/7/1  1999/8/1  1999/9/1  1999/10/1
## [57] 1999/11/1 1999/12/1 2000/1/1 
## 59 Levels: 1995/10/1 1995/11/1 1995/12/1 1995/3/1 1995/4/1 ... 2000/1/1
library(xts)
library(PerformanceAnalytics)
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
date1 = as.Date(firm_data1[,1], "%Y/%m/%d")
#convert firm_data1 into time series data: xts
firm_data1.xts = as.xts(firm_data1[,-1], order.by = date1)
n <- dim(firm_data1.xts)[2]
constraints = new.constraints(n, lb = -Inf, ub = +Inf)
# SUM x.i = 1
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)  
ia <- create.historical.ia(firm_data1.xts, 12)
names(ia)
##  [1] "hist.returns"      "nperiod"           "index"            
##  [4] "n"                 "symbols"           "risk"             
##  [7] "correlation"       "cov"               "expected.return"  
## [10] "annual.factor"     "arithmetic.return" "geometric.return"
#s0 <- apply(coredata(firm_data1.xts),2,sd)     
#ia$cov <- cor(coredata(firm_data1.xts), use='complete.obs',method='pearson') * (s0 %*% t(s0))
weight <- min.risk.portfolio(ia, constraints)
## Loading required package: kernlab
## 
## Attaching package: 'kernlab'
## The following object is masked _by_ '.GlobalEnv':
## 
##     cross
weight
## [1] 0.3635998 0.1936537 0.4427465
#===================================================
# Try to plot efficient frontier using SIT
# create sample historical input assumptions
ia <- create.historical.ia(firm_data1.xts, 12) # 12 is annual factor for monthly return
# create long-only, fully invested efficient frontier
# 0 <= x.i <= 1
# If short sale allowed: constraints = new.constraints(n, lb = -Inf, ub = +Inf)
constraints = new.constraints(n, lb = 0, ub = 1)
constraints = add.constraints(diag(n), type='>=', b=0, constraints)
constraints = add.constraints(diag(n), type='<=', b=1, constraints)
# SUM x.i = 1
constraints = add.constraints(rep(1, n), 1, type = '=', constraints)
#
weight <- min.risk.portfolio(ia, constraints)
# create efficient frontier
ifelse(!require(corpcor), install.packages("corpcor"), library(corpcor))
## Loading required package: corpcor
## 
## Attaching package: 'corpcor'
## The following object is masked _by_ '.GlobalEnv':
## 
##     cov.shrink
## [1] "corpcor"
ifelse(!require(lpSolve), install.packages("lpSolve"), library(lpSolve))
## Loading required package: lpSolve
## [1] "lpSolve"
ef = portopt(ia, constraints, 50, 'Efficient Frontier') 
ef
## $weight
##         Nordstrom  Starbucks Microsoft
##  [1,] 0.363599832 0.19365370 0.4427465
##  [2,] 0.354174974 0.19474884 0.4510762
##  [3,] 0.344750115 0.19584397 0.4594059
##  [4,] 0.335325256 0.19693910 0.4677356
##  [5,] 0.325900397 0.19803423 0.4760654
##  [6,] 0.316475538 0.19912936 0.4843951
##  [7,] 0.307050679 0.20022449 0.4927248
##  [8,] 0.297625821 0.20131963 0.5010546
##  [9,] 0.288200962 0.20241476 0.5093843
## [10,] 0.278776103 0.20350989 0.5177140
## [11,] 0.269351244 0.20460502 0.5260437
## [12,] 0.259926385 0.20570015 0.5343735
## [13,] 0.250501527 0.20679529 0.5427032
## [14,] 0.241076668 0.20789042 0.5510329
## [15,] 0.231651809 0.20898555 0.5593626
## [16,] 0.222226950 0.21008068 0.5676924
## [17,] 0.212802091 0.21117581 0.5760221
## [18,] 0.203377232 0.21227094 0.5843518
## [19,] 0.193952374 0.21336608 0.5926816
## [20,] 0.184527515 0.21446121 0.6010113
## [21,] 0.175102656 0.21555634 0.6093410
## [22,] 0.165677797 0.21665147 0.6176707
## [23,] 0.156252938 0.21774660 0.6260005
## [24,] 0.146828079 0.21884173 0.6343302
## [25,] 0.137403221 0.21993687 0.6426599
## [26,] 0.127978362 0.22103200 0.6509896
## [27,] 0.118553503 0.22212713 0.6593194
## [28,] 0.109128644 0.22322226 0.6676491
## [29,] 0.099703785 0.22431739 0.6759788
## [30,] 0.090278927 0.22541252 0.6843085
## [31,] 0.080854068 0.22650766 0.6926383
## [32,] 0.071429209 0.22760279 0.7009680
## [33,] 0.062004350 0.22869792 0.7092977
## [34,] 0.052579491 0.22979305 0.7176275
## [35,] 0.043154632 0.23088818 0.7259572
## [36,] 0.033729774 0.23198331 0.7342869
## [37,] 0.024304915 0.23307845 0.7426166
## [38,] 0.014880056 0.23417358 0.7509464
## [39,] 0.005455197 0.23526871 0.7592761
## [40,] 0.000000000 0.22636817 0.7736318
## [41,] 0.000000000 0.20373135 0.7962686
## [42,] 0.000000000 0.18109453 0.8189055
## [43,] 0.000000000 0.15845772 0.8415423
## [44,] 0.000000000 0.13582090 0.8641791
## [45,] 0.000000000 0.11318408 0.8868159
## [46,] 0.000000000 0.09054727 0.9094527
## [47,] 0.000000000 0.06791045 0.9320895
## [48,] 0.000000000 0.04527363 0.9547264
## [49,] 0.000000000 0.02263682 0.9773632
## [50,] 0.000000000 0.00000000 1.0000000
## 
## $return
##            [,1]
##  [1,] 0.3729416
##  [2,] 0.3786335
##  [3,] 0.3843255
##  [4,] 0.3900175
##  [5,] 0.3957094
##  [6,] 0.4014014
##  [7,] 0.4070934
##  [8,] 0.4127853
##  [9,] 0.4184773
## [10,] 0.4241693
## [11,] 0.4298613
## [12,] 0.4355532
## [13,] 0.4412452
## [14,] 0.4469372
## [15,] 0.4526291
## [16,] 0.4583211
## [17,] 0.4640131
## [18,] 0.4697050
## [19,] 0.4753970
## [20,] 0.4810890
## [21,] 0.4867809
## [22,] 0.4924729
## [23,] 0.4981649
## [24,] 0.5038568
## [25,] 0.5095488
## [26,] 0.5152408
## [27,] 0.5209327
## [28,] 0.5266247
## [29,] 0.5323167
## [30,] 0.5380087
## [31,] 0.5437006
## [32,] 0.5493926
## [33,] 0.5550846
## [34,] 0.5607765
## [35,] 0.5664685
## [36,] 0.5721605
## [37,] 0.5778524
## [38,] 0.5835444
## [39,] 0.5892364
## [40,] 0.5949283
## [41,] 0.6006203
## [42,] 0.6063123
## [43,] 0.6120042
## [44,] 0.6176962
## [45,] 0.6233882
## [46,] 0.6290801
## [47,] 0.6347721
## [48,] 0.6404641
## [49,] 0.6461561
## [50,] 0.6518480
## 
## $risk
##  [1] 0.2539256 0.2539584 0.2540569 0.2542210 0.2544506 0.2547454 0.2551053
##  [8] 0.2555300 0.2560192 0.2565724 0.2571893 0.2578695 0.2586123 0.2594174
## [15] 0.2602840 0.2612117 0.2621997 0.2632473 0.2643540 0.2655188 0.2667411
## [22] 0.2680201 0.2693549 0.2707448 0.2721888 0.2736862 0.2752361 0.2768376
## [29] 0.2784897 0.2801917 0.2819425 0.2837413 0.2855873 0.2874794 0.2894167
## [36] 0.2913985 0.2934237 0.2954915 0.2976009 0.2997974 0.3024626 0.3056698
## [43] 0.3094021 0.3136409 0.3183658 0.3235556 0.3291882 0.3352415 0.3416929
## [50] 0.3485204
## 
## $name
## [1] "Efficient Frontier"
plot.ef(ia, list(ef), transition.map=F)  

# find maximum sharpe portfolio
max(portfolio.return(ef$weight,ia) /  portfolio.risk(ef$weight,ia))
## [1] 1.985767
# plot minimum variance portfolio
weight = min.var.portfolio(ia,constraints)  
points(100 * portfolio.risk(weight,ia), 100 * portfolio.return(weight,ia), pch=15, col='red')
portfolio.return(weight,ia) /  portfolio.risk(weight,ia)
##          [,1]
## [1,] 1.468704
# plot maximum Sharpe or tangency portfolio
weight = max.sharpe.portfolio()(ia,constraints) 
points(100 * portfolio.risk(weight,ia), 100 * portfolio.return(weight,ia), pch=15, col='orange')
portfolio.return(weight,ia) /  portfolio.risk(weight,ia)
##          [,1]
## [1,] 1.985799
plota.legend('Minimum Variance,Maximum Sharpe','red,orange', x='topright')