SCM

SCM Repository

[rmetrics] Annotation of /pkg/fPortfolio/R/102C-MarkowitzPortfolio.R
ViewVC logotype

Annotation of /pkg/fPortfolio/R/102C-MarkowitzPortfolio.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 61 - (view) (download)

1 : wuertz 17
2 :     # This library is free software; you can redistribute it and/or
3 :     # modify it under the terms of the GNU Library General Public
4 :     # License as published by the Free Software Foundation; either
5 :     # version 2 of the License, or (at your option) any later version.
6 :     #
7 :     # This library is distributed in the hope that it will be useful,
8 :     # but WITHOUT ANY WARRANTY; without even the implied warranty of
9 :     # MERCHANTABILITY or FITNESS FOR A PARTICULAR Description. See the
10 :     # GNU Library General Public License for more details.
11 :     #
12 :     # You should have received a copy of the GNU Library General
13 :     # Public License along with this library; if not, write to the
14 :     # Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
15 :     # MA 02111-1307 USA
16 :    
17 :     # Copyrights (C)
18 :     # for this R-port:
19 :     # 1999 - 2004, Diethelm Wuertz, GPL
20 :     # Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
21 :     # info@rmetrics.org
22 :     # www.rmetrics.org
23 :     # for the code accessed (or partly included) from other R-ports:
24 :     # see R's copyright and license files
25 :     # for the code accessed (or partly included) from contributed R-ports
26 :     # and other sources
27 :     # see Rmetrics's copyright file
28 :    
29 :    
30 :     ################################################################################
31 :     # FUNCTION: EFFICIENT FRONTIER:
32 :     # fPFOLIO S4 Portfolio Class
33 :     # frontierMarkowitz Efficient frontier of mean-var portfolio
34 :     # montecarloMarkowitz Adds randomly created portfolios
35 :     # portfolioMarkowitz Mean-variance Markowitz (target) portfolio
36 :     # METHODS: DESCRIPTION:
37 : wuertz 52 # .frontier.default S3: Extract points on the efficient frontier
38 : wuertz 17 # print.fPFOLIO S3: Print method for objects of class fPFOLIO
39 :     # plot.fPFOLIO S3: Plot method for objects of class fPFOLIO
40 :     # summary.fPFOLIO S3: Summary method for objects of class fPFOLIO
41 :     # FUNCTIONS: SPECIAL PORTFOLIOS:
42 :     # .tangencyMarkowitz Adds tangency portfolio
43 :     # .equalweightsMarkowitz Adds equal weights Portfolio
44 :     # .frontierShortSelling... Efficient Frontier of Short Selling Markowitz
45 :     # FUNCTIONS: BUILTIN FROM PACKAGE:
46 :     # .portfolio.optim Function from R-package tseries
47 :     # .solve.QP Function from R-package quadprog
48 :     ################################################################################
49 :    
50 :    
51 :    
52 :     setClass("fPFOLIO",
53 :     representation(
54 :     call = "call",
55 :     method = "character",
56 :     model = "character",
57 :     data = "list",
58 :     pfolio = "list",
59 :     title = "character",
60 :     description = "character")
61 :     )
62 :    
63 :    
64 :     # ------------------------------------------------------------------------------
65 :    
66 :    
67 :     frontierMarkowitz =
68 :     function(data, Rf = 0, short = FALSE, length = 300, r.range = NULL,
69 :     s.range = NULL, title = NULL, description = NULL, ...)
70 :     { # A function implemented by Diethelm Wuertz
71 :    
72 :     # Description:
73 :     # Calculates the efficient frontier from a matrix
74 :     # of either market or simulated assets given in matrix "x".
75 :     # Each time series represents a column in this matrix.
76 :    
77 :     # Arguments:
78 :     # data - either a rectangular time series object which can be
79 :     # transformed into a matrix or a list with two named entries,
80 :     # mu - the returns of the multivariate series, and
81 :     # Sigma - the covariance matrix of the multivariate series.
82 :     # Rf - risk free rate, a numerical value by default 0.
83 :     # short - a logical flag, if set to TRUE short selling would be
84 :     # allowed otherwise not wich is the default setting.
85 :     # length - Number of equidistant return points on the Efficient
86 :     # frontier.
87 :     # r.range - Plot range of returns, if NULL, the range will be created
88 :     # automatically.
89 :     # s.range - Plot range of standard deviations, if NULL, the range
90 :     # will be created automatically.
91 :    
92 :     # Value:
93 :     # data - the matrix of asset time series
94 :     # ps - the portfolio standard deviation on the efficient frontier
95 :     # pm - the portfolio means on the efficient frontier
96 :     # cov - the covariance matrix of the assets
97 :     # r.range - the plot range for the mean returns, y-axis
98 :     # s.range - the plot range for the standard deviations, x-Axis
99 :     # diversification - number of portfolio weights on the efficient frontier
100 :     # call - the call string
101 :     # type - the type of portfolio optimizatiom applied
102 :     # series - the name of asset time series
103 :     # Rf, Rm, Sm, t.weights - not used, added by tangencyPortfolio
104 :    
105 :     # Notes:
106 :     # * Portfolio optimization is done by a method specified by the
107 :     # argument type, at the moment this is "markowitz". No other
108 :     # methods are available so far.
109 :     # * The plot range is determined by default automatically from
110 :     # the market time series through r,range and s.range, if not
111 :     # otherwise specified.
112 :     # * The function returns an object of class portfolio with several
113 :     # components stored in a list.
114 :    
115 :     # FUNCTION:
116 :    
117 :     # Short Selling Allowed ?
118 :     if (short) {
119 : wuertz 52 return(.frontierShortSellingMarkowitz(data = data, Rf = Rf,
120 :     length = length, r.range = r.range, s.range = s.range,
121 :     title = title, description = description, ...))
122 : wuertz 17 }
123 :    
124 :     # Check Data:
125 :     if (is.list(data)) {
126 : wuertz 49 # data is a list with entries list(mu, Sigma)
127 :     x.mat = NA
128 : wuertz 17 mu = data$mu
129 :     Sigma = data$Sigma
130 :     plottype = "mucov"
131 :     } else {
132 : wuertz 49 # data is a rectangular time series object:
133 : wuertz 17 x.mat = as.matrix(data)
134 :     # Mean Vector and Covariance Matrix:
135 :     mu = apply(x.mat, MARGIN = 2, FUN = mean)
136 :     Sigma = apply(x.mat, MARGIN = 2, FUN = var)
137 :     plottype = "series"
138 :     }
139 :    
140 :     # Ranges for mean and Standard Deviation:
141 :     if (is.null(r.range)) r.range = range(mu)
142 :     if (is.null(s.range)) s.range = c(0, max(sqrt(diag(Sigma))))
143 :    
144 :     # Calculate Efficient Frontier:
145 :     pmmin = min(r.range)
146 :     pmmax = max(r.range)
147 : wuertz 49 ps = pm = diversification = rep(0, times = length)
148 : wuertz 17 eps = 1.0e-6
149 :     k = 0
150 : wuertz 49 for (PM in seq(pmmin+eps, pmmax-eps, length = length)) {
151 : wuertz 17 k = k+1
152 : wuertz 49 ef = .portfolio.optim(pm = PM, returns = mu, covmat = Sigma)
153 : wuertz 17 pm[k] = ef$pm
154 :     ps[k] = ef$ps
155 :     diversification[k] = length (ef$pw[ef$pw > 0.01])
156 :     }
157 :    
158 :     # Result:
159 :     pfolio = list(
160 :     what = "frontier",
161 :     plottype = plottype,
162 : wuertz 49 data = data,
163 : wuertz 17 pw = NA, pm = pm, ps = ps,
164 : wuertz 49 returns = mu, cov = Sigma,
165 : wuertz 17 r.range = r.range, s.range = s.range,
166 :     Rf = NA, Rm = NA, Sm = NA,
167 :     Rew = NA, Sew = NA,
168 :     t.weights = NA,
169 :     diversification = diversification)
170 :    
171 :     # Add Tangency Portfolio:
172 :     pfolio = .tangencyMarkowitz(pfolio, Rf = Rf, add = FALSE)
173 :    
174 :     # Add Equal Weights Portfolio:
175 :     pfolio = .equalweightsMarkowitz(pfolio, add = FALSE)
176 :    
177 :     # Title:
178 :     if (is.null(title))
179 :     title = "Mean-Variance Portfolio Optimization"
180 :    
181 :     # Description:
182 :     if (is.null(description))
183 :     description = as.character(date())
184 :    
185 :     # Return Value:
186 :     new ("fPFOLIO",
187 :     call = as.call(match.call()),
188 :     method = "Quadratic Programming",
189 :     model = "Markowitz Portfolio",
190 :     data = list(x = data),
191 :     pfolio = pfolio,
192 :     title = as.character(title),
193 :     description = as.character(description)
194 :     )
195 :     }
196 :    
197 :    
198 :     # ------------------------------------------------------------------------------
199 :    
200 :    
201 :     montecarloMarkowitz =
202 :     function(object, mc = 5000, doplot = FALSE, add = TRUE, ...)
203 :     { # A function implemented by Diethelm Wuertz
204 :    
205 :     # Description:
206 :     # This function adds by Monte carlo simulation to an effiicient
207 :     # frontier plot randomly selected portfolios and determines and
208 :     # plots a MC efficient frontier.
209 :    
210 :     # Arguments:
211 :     # object - An object of class portfolio
212 :     # mc - Number of Monte Carlo steps
213 :    
214 :     # Value:
215 :     # An updated object of class portfolio with the following
216 :     # components added/updated:
217 :     # mcpm - Monte Carlo mean returns on MC efficient frontier
218 :     # mcps - Monte Carlo standard deviations on MC efficient forntier
219 :    
220 :     # FUNCTION:
221 :    
222 :     # Extract Portfolio:
223 :     keep = object
224 :     object = object@pfolio
225 :    
226 :     # Internal function:
227 :     records =
228 :     function (data, conf.level = 0.95, ...) {
229 :     data = as.numeric(data)
230 :     record = cummax(data)
231 :     expected = cumsum(1/(1:length(data)))
232 :     se = sqrt(expected - cumsum(1/((1:length(data))^2)))
233 :     trial = (1:length(data))[!duplicated(record)]
234 :     record = unique(record)
235 :     number = 1:length(record)
236 :     expected = expected[trial]
237 :     se = se[trial]
238 :     data.frame(number, record, trial, expected, se)
239 :     }
240 :    
241 :     # Settings:
242 :     returns = object$returns
243 :     covar = object$cov
244 :     dimension = length(returns)
245 :     r.range = object$r.range
246 :     s.range = object$s.range
247 :    
248 :     # Plot frame:
249 :     if (doplot) {
250 :     plot(x = s.range, y = r.range, type = "n",
251 :     xlab = "Standard Deviation", ylab = "Return", ...)
252 :     }
253 :    
254 :     # Monte Carlo Loop:
255 :     rr = ss = rep(0, mc)
256 :     for (k in 1:mc) {
257 :     weights = abs(rcauchy(dimension))
258 :     weights = weights/sum(weights)
259 :     rp = weights %*% returns
260 :     sp = weights %*% covar %*% weights
261 :     sp = sqrt(sp)
262 :     if (add) points(sp, rp, cex = 0.1)
263 :     ss[k] = sp
264 :     rr[k] = rp
265 :     }
266 :    
267 :     # Simulated Efficient Frontier:
268 :     # Lower Part:
269 :     rr = -rr
270 :     order = order(ss)
271 :     xsp.lower = ss[order]
272 :     xrp.lower = rr[order]
273 :     trial = records(xrp.lower)$trial
274 :     lines(xsp.lower[trial], -xrp.lower[trial], col = "red")
275 :     # Upper Part:
276 :     rr = -rr
277 :     order = order(ss)
278 :     xsp.upper = ss[order]
279 :     xrp.upper = rr[order]
280 :     trial = records(xrp.upper)$trial
281 :     lines(xsp.upper[trial], xrp.upper[trial], col = "red")
282 :     lines(x = s.range, y = mean(rr)*c(1,1), col = "green")
283 :    
284 :     # Renew Asset Texts:
285 :     if (add) {
286 :     points(sqrt(diag(covar)), returns, col = "red")
287 :     s.delta = diff(s.range)/50
288 :     text(sqrt(diag(covar))+s.delta, returns, as.character(1:dimension),
289 :     col = "red")
290 :     }
291 :    
292 :     # Return Monte Carlo Frontier:
293 :     object$mcpm = c(rev(xrp.lower[trial]), xrp.upper[trial])
294 :     object$mcps = c(rev(xsp.lower[trial]), xrp.upper[trial])
295 :    
296 :     # Return Value
297 :     invisible(keep)
298 :     }
299 :    
300 :    
301 :     # ------------------------------------------------------------------------------
302 :    
303 :    
304 :     portfolioMarkowitz =
305 :     function(x, targetReturn, title = NULL, description = NULL)
306 :     { # A function implemented by Diethelm Wuertz
307 :    
308 :     # Description:
309 :     # Optimizes a mean-var portfolio for a given desired return
310 :    
311 :     # Arguments:
312 :     # x - portfolio of assets, a matrix or an object which can be
313 :     # transformed to a matrix
314 :     # targetReturn - the desired return, pm must be smaller than the
315 :     # maximum and larger than the minimum asset return. Short
316 :     # selling is not allowed.
317 :    
318 :     # FUNCTION:
319 :    
320 :     # Transform to matrix:
321 :     x = as.matrix(x)
322 : wuertz 49 MU = apply(x, 2, mean)
323 : wuertz 17 COV = cov(x)
324 :    
325 :     # Quadratic Programming:
326 :     pfolio = list()
327 :     pfolio$what = "portfolio"
328 :     pfolio$method = "QP"
329 : wuertz 49 opt = .portfolio.optim(pm = targetReturn, returns = MU, covmat = COV)
330 : wuertz 17 pfolio$pw = opt$pw
331 :     pfolio$pm = opt$pm
332 :     pfolio$ps = opt$ps
333 :     pfolio$returns = apply(x, 2, mean)
334 :     pfolio$cov = COV
335 :     pfolio$r.range = range(apply(x, 2, mean))
336 :     pfolio$s.range = c(0, max(sqrt(diag(COV))))
337 :    
338 :     # Title:
339 :     if (is.null(title))
340 :     title = "Mean-Variance Portfolio Optimization"
341 :    
342 :     # Description:
343 :     if (is.null(description))
344 :     description = as.character(date())
345 :    
346 :     # Return Value:
347 :     new ("fPFOLIO",
348 :     call = as.call(match.call()),
349 :     method = "Quadratic Programming",
350 :     model = "Markowitz Portfolio",
351 :     data = list(x = x),
352 :     pfolio = pfolio,
353 :     title = as.character(title),
354 :     description = as.character(description)
355 :     )
356 :     }
357 :    
358 :    
359 :     # ******************************************************************************
360 :    
361 :    
362 : wuertz 52 .frontier =
363 :     function (object, ...)
364 :     { # A function implemented by Diethelm Wuertz
365 :    
366 :     # FUNCTION:
367 :    
368 :     # Return Value:
369 :     UseMethod(".frontier")
370 :     }
371 :    
372 :    
373 :     # ------------------------------------------------------------------------------
374 :    
375 :    
376 :     .frontier.default =
377 : wuertz 61 function(object, ...)
378 : wuertz 17 { # A function implemented by Diethelm Wuertz
379 :    
380 :     # Description:
381 :     # S3 method for an object of class "fPFOLIO" to extract
382 :     # the points on the efficient frontier
383 :    
384 :     # FUNCTION:
385 :    
386 :     # Frontier:
387 :     pfolio = object@pfolio
388 :     mu = pfolio$pm
389 : wuertz 52 sigma = pfolio$ps
390 : wuertz 17
391 :     # Return Value:
392 : wuertz 52 data.frame(sigma = sigma, mu = mu)
393 : wuertz 17 }
394 :    
395 :    
396 :     # ******************************************************************************
397 :    
398 :    
399 :     print.fPFOLIO =
400 :     function(x, ...)
401 :     { # A function implemented by Diethelm Wuertz
402 :    
403 :     # Description:
404 :     # S3 Print Method for an object of class "fPFOLIO"
405 :    
406 :     # Arguments:
407 :     # x - an object of class "fPFOLIO"
408 :    
409 :     # FUNCTION:
410 :    
411 :     # Extract Portfolio:
412 :     pfolio = x@pfolio
413 :    
414 :     # Title:
415 :     cat("\nTitle:\n ")
416 :     cat(as.character(x@title), "\n")
417 :    
418 :     # Call:
419 :     cat("\nCall:\n ")
420 :     print.default(x@call)
421 :    
422 :     # Portfolio Optimization:
423 :     if (pfolio$what == "portfolio") {
424 :     # Weights:
425 :     weights = round(pfolio$pw, digits = 2)
426 :     names(weights) = as.character(1:length(weights))
427 :     cat("\nPortfolio Weights:\n ")
428 :     print(weights[weights > 0])
429 :     # Check Sum of Weights:
430 :     cat("\nSum of Weights:\n ")
431 :     print(round(sum(pfolio$pw), digits = 2))
432 :     # Target Returns:
433 :     cat("\nTarget Return(s):\n ")
434 :     print(round(pfolio$pm, digits = 4))
435 :     # Target Risk:
436 :     cat("\nTarget Risk(s):\n ")
437 :     print(round(pfolio$ps, digits = 4))
438 :     }
439 :    
440 :     # Frontier:
441 :     if (pfolio$what == "frontier") {
442 :     # Efficient Frontier:
443 :     cat("\nEfficient Frontier - Returns:\n ")
444 :     r.digits = abs(round(log10(max(pfolio$pm)), 0)) + 4
445 :     print(round(pfolio$pm, digits = r.digits))
446 :     cat("\nEfficient Frontier - Standard Deviations:\n ")
447 :     s.digits = abs(round(log10(max(pfolio$ps)), 0)) + 4
448 :     print(round(pfolio$ps, digits = s.digits))
449 :     # Tangency Portfolio:
450 :     cat("\nTangency Portfolio:\n ")
451 :     cat("Risk free rate:", round(pfolio$Rf, digits = r.digits),
452 :     " Return:", round(pfolio$Rm, digits = r.digits),
453 :     " Risk:", round(pfolio$Sm, digits = s.digits), "\n")
454 :     # Equal Weights Portfolio:
455 :     cat("\nEqual Weights Portfolio:\n ")
456 :     cat("Risk:", round(pfolio$Rew, digits = r.digits),
457 :     " Return:", round(pfolio$Sew, digits = s.digits), "\n")
458 :     }
459 :    
460 :     # Description:
461 :     cat("\nDescription:\n ")
462 :     print(x@description)
463 :    
464 :     # Return Value:
465 :     invisible(x)
466 :     }
467 :    
468 :    
469 :     # ------------------------------------------------------------------------------
470 :    
471 :    
472 :     plot.fPFOLIO =
473 :     function(x, alpha = 0.05, mc = 500, which = "ask", ...)
474 :     { # A function implemented by Diethelm Wuertz
475 :    
476 :     # Description:
477 :     # S3 Print Method for an object of class "fPFOLIO"
478 :    
479 :     # Arguments:
480 :     # x - an object of class "fPFOLIO"
481 :    
482 :     # FUNCTION:
483 :    
484 :     # Check:
485 :     if (as.character(x@call[1]) == "portfolioMarkowitz")
486 :     stop("Objects returned by portfolioMarkowitz() cannot be plotted")
487 :    
488 :     # Settings:
489 :     .mcSteps <<- mc
490 :     .alphaLevel <<- alpha
491 :     if (which[1] == "ask") {
492 :     .interactive <<- TRUE
493 :     } else {
494 :     .interactive <<- FALSE
495 :     }
496 :    
497 :     # Plot Efficient Frontier:
498 :     plot.1 <<- function(x) {
499 :     pfolio = x@pfolio
500 :     # Extract:
501 :     dim = length(pfolio$returns)
502 :     returns = pfolio$returns
503 :     covar = pfolio$cov
504 :     r.range = pfolio$r.range
505 :     s.range = pfolio$s.range
506 :     diversification = pfolio$diversification
507 :     # 1. Plot Graph Frame:
508 :     plot(x = s.range, y = r.range, type = "n",
509 :     xlab = "Standard Deviation", ylab = "Return",
510 :     main = "Mean Variance Portfolio")
511 :     # 2. Plot Individual Assets:
512 :     points(sqrt(diag(covar)), returns, pch = 20, col = "red")
513 :     risk = sqrt(diag(covar)) + diff(s.range) / 50
514 :     text(risk, returns, as.character(1:dim), col = "red")
515 :     # 3. Plot Efficient Frontier:
516 :     points(pfolio$ps, pfolio$pm, col = "blue", cex = 0.1)
517 :     # 4. Add Tangency Portfolio:
518 :     # Find the Index (location of the Tangency Point):
519 :     delta.pm = diff(pfolio$pm)
520 :     delta.ps = diff(pfolio$ps)
521 :     slope1 = slope.keep = delta.pm/delta.ps
522 :     slope1 = slope1[slope.keep > 0]
523 :     slope2 = (pfolio$pm-pfolio$Rf) / pfolio$ps
524 :     diff.slope2 = diff(slope2) / 2
525 :     slope2 = slope2[1:length(diff.slope2)] + diff.slope2
526 :     slope2 = slope2[slope.keep > 0]
527 :     cross = abs(slope2 - slope1)
528 :     index = order(cross)[1]
529 :     pm.pos = pfolio$pm[slope.keep > 0]
530 :     ps.pos = pfolio$ps[slope.keep > 0]
531 :     # Here are the mean return and the standard deviation
532 :     mean.m = (pm.pos[index] + pm.pos[index+1]) / 2
533 :     sigma.m = (ps.pos[index] + ps.pos[index+1]) / 2
534 :     # ... together with the slope
535 :     slope = (mean.m - pfolio$Rf) / sigma.m
536 :     # Add the Capital Market Line and the Tangency Point to the plot
537 :     lines(c(0, max(pfolio$ps)), c(pfolio$Rf, slope* max(pfolio$ps) +
538 :     pfolio$Rf), lwd = 2 )
539 :     points(sigma.m, slope*sigma.m + pfolio$Rf, col = "green",
540 :     pch = 20, cex = 1.5)
541 :     s.delta = max(pfolio$ps) / 10
542 :     text(sigma.m - s.delta, slope*sigma.m + pfolio$Rf, "CML")
543 :     # 5. Add Equal Weights Portfolio:
544 :     # Calculate equally weighted portfolio:
545 :     ew.weights = rep(1/dim, times = dim)
546 :     Rew = (ew.weights %*% returns)
547 :     Sew = (ew.weights %*% covar %*% ew.weights)
548 :     points(sqrt(Sew), Rew, pch = 20, col = "steelblue4", cex = 1.5)
549 :     }
550 :     # ... add Monte Carlo Portfolios:
551 :     plot.2 <<- function(x) {
552 :     plot.1(x)
553 :     if (.interactive) {
554 :     mc = readline("Number of Monte Carlo Steps > ")
555 :     if (mc != "") .mcSteps <<- as.integer(mc)
556 :     cat("\nTotal Number of Monte Carlo Steps:", .mcSteps, "\n")
557 :     }
558 :     tmp = montecarloMarkowitz(x, mc = .mcSteps, doplot = FALSE, add = TRUE)
559 :    
560 :     }
561 :     # Bar Plot of Tangency Weights:
562 :     plot.3 <<- function(x) {
563 :     pfolio = x@pfolio
564 :     weights = pfolio$t.weights ## [pfolio$t.weights > 0]
565 :     ## pie(weights, col = rainbow(length(weights)), radius = 0.9,
566 :     ## main = "Tangency Portfolio Weights")
567 :     barplot(weights, cex.names = 0.7, col = "steelblue")
568 :     abline (h = 0)
569 :     }
570 :     # Plot Diversification:
571 :     plot.4 <<- function(x) {
572 :     pfolio = x@pfolio
573 :     if (is.na(pfolio$diversification)) {
574 :     cat ("\nNo time series data available.\n")
575 :     } else {
576 :     xx = pfolio$pm
577 :     yy = pfolio$diversification
578 :     data = data.frame(cbind(xx, yy))
579 :     plot(pfolio$pm, pfolio$diversification, type = "l", xlab = "Return",
580 :     ylab = "Number of Assets", main = "Diversification")
581 :     zz = loess(yy ~ xx, data)
582 :     lines(xx, zz$fit, col = "blue")
583 :     points(pfolio$Rm, length(weights), col = "green", cex = 1.5)
584 :     }
585 :     }
586 :     # 5. Plot Cumulated Return:
587 :     plot.5 <<- function(x) {
588 :     pfolio = x@pfolio
589 :     if (pfolio$plottype == "mucov") {
590 :     cat ("\nNo time series data available.\n")
591 :     } else {
592 :     plot(cumsum(as.matrix(pfolio$data) %*% pfolio$t.weights),
593 :     type = "l", col = "steelblue",
594 :     ylab = "Cumulated Return", main = "Portfolio Return")
595 :     grid()
596 :     }
597 :     }
598 :     # 6. Portfolio Risk Histogram:
599 :     plot.6 <<- function(x) {
600 :     pfolio = x@pfolio
601 :     if (pfolio$plottype == "mucov") {
602 :     cat ("\nNo time series data available.\n")
603 :     } else {
604 :     weights = pfolio$t.weights[pfolio$t.weights > 0]
605 :     pfolioHist(x@data, weights = weights, alpha = .alphaLevel)
606 :     }
607 :     }
608 :    
609 :     # Plot:
610 :     interactivePlot(
611 :     x = x,
612 :     choices = c(
613 :     "Efficient Frontier",
614 :     "... with Monte Carlo Portfolios",
615 :     "Bar Chart of Tangency Portfolio",
616 :     "Diversification Plot",
617 :     "Cumulated Return Tangency Portfolio",
618 :     "Risk of Tangency Portfolio"),
619 :     plotFUN = c(
620 :     "plot.1",
621 :     "plot.2",
622 :     "plot.3",
623 :     "plot.4",
624 :     "plot.5",
625 :     "plot.6"),
626 :     which = which)
627 :    
628 :    
629 :     # Return Value:
630 :     invisible(x)
631 :     }
632 :    
633 :    
634 :     # ------------------------------------------------------------------------------
635 :    
636 :    
637 :     summary.fPFOLIO =
638 :     function(object, ...)
639 :     { # A function implemented by Diethelm Wuertz
640 :    
641 :     # Description:
642 :     # Summary Method for an object of class "fPFOLIO"
643 :    
644 :     # FUNCTION:
645 :    
646 :     # Extract Portfolio:
647 :     pfolio = object@pfolio
648 :    
649 :     # Print Result:
650 :     print(object)
651 :    
652 :     # Plot:
653 :     plot(object, ...)
654 :    
655 :     # Return Value:
656 :     invisible(object)
657 :     }
658 :    
659 :    
660 :     ################################################################################
661 :    
662 :    
663 :     .tangencyMarkowitz =
664 :     function(object, Rf = 0, add = TRUE)
665 :     { # A function implemented by Diethelm Wuertz
666 :    
667 :     # Description:
668 :     # This function calculates the properties of the tangency portfolio,
669 :     # and adds the "capital market line" and the tangency portfolio, Rm
670 :     # and Sm, to the efficient frontier.
671 :    
672 :     # NOTE:
673 :     # Not for all risk free rates there exists a tangency portfolio!!!
674 :    
675 :     # Arguments:
676 :     # object - An object of class portfolio
677 :    
678 :     # Value:
679 :     # An updated object of class portfolio with the following
680 :     # components added/updated:
681 :     # Rf - risk free rate of return
682 :     # Rm - rate of return of the tangency portfolio
683 :     # Sm - standardard deviation of the tangency portfolio
684 :     # t.weights - weights of the tangency portfolio
685 :    
686 :     # Note:
687 :     # I have not yet checked if this works for all values of Rf !!
688 :     # I hope it works properly even if there exists no tangency portfolio.
689 :    
690 :     # FUNCTION:
691 :    
692 :     # Find the Index (location of the Tangency Point:
693 :     delta.pm = diff(object$pm)
694 :     delta.ps = diff(object$ps)
695 :     slope1 = slope.keep = delta.pm/delta.ps
696 :     slope1 = slope1[slope.keep > 0]
697 :     slope2 = (object$pm-Rf) / object$ps
698 :     diff.slope2 = diff(slope2) / 2
699 :     slope2 = slope2[1:length(diff.slope2)] + diff.slope2
700 :     slope2 = slope2[slope.keep > 0]
701 :     cross = abs(slope2 - slope1)
702 :     index = order(cross)[1]
703 :     pm.pos = object$pm[slope.keep > 0]
704 :     ps.pos = object$ps[slope.keep > 0]
705 :    
706 :     # Here are the mean return and the standard deviation
707 :     mean.m = (pm.pos[index] + pm.pos[index+1]) / 2
708 :     sigma.m = (ps.pos[index] + ps.pos[index+1]) / 2
709 :     # ... together with the slope
710 :     slope = (mean.m - Rf) / sigma.m
711 :     x1 = 0
712 :     y1 = Rf
713 :     x2 = max(object$ps)
714 :     y2 = slope*x2 + Rf
715 :     xm = sigma.m
716 :     ym = slope*sigma.m + Rf
717 :    
718 :     # Add the Capital Market Line and the Tangency Point to the plot
719 :     if (add) {
720 :     lines(c(x1, x2), c(y1, y2), lwd = 2 )
721 :     points(xm, ym, col = "green", cex = 1.5)
722 :     s.delta = max(object$ps) / 10
723 :     text(xm - s.delta, ym, "CML")}
724 :    
725 :     # Calculate Tangency Portfolio, if it exists!
726 :     # The test is true, when the slope could be evaluated.
727 :     if (!is.na(slope)) {
728 : wuertz 49 tangency = .portfolio.optim(pm = ym, returns = object$returns,
729 : wuertz 17 covmat = object$cov)
730 :     # Portfolio Weights in Percent:
731 :     pw = 100 * tangency$pw
732 :     # Remove those smaller 1 %% Promille!!
733 :     pw = round(pw*(sign(pw - 0.1) + 1)/2, digits = 2)
734 :     object$t.weights = pw/100
735 :     }
736 :    
737 :     # Add to Output:
738 :     object$Rf = Rf
739 :     object$Rm = ym
740 :     object$Sm = xm
741 :    
742 :     # Return Value:
743 :     object
744 :     }
745 :    
746 :    
747 :     # ------------------------------------------------------------------------------
748 :    
749 :    
750 :     .equalweightsMarkowitz =
751 :     function(object, add = TRUE)
752 :     { # A function implemented by Diethelm Wuertz
753 :    
754 :     # Description:
755 :     # This function adds the equal weighted Portfolio to
756 :     # efficient frontier plot.
757 :    
758 :     # FUNCTION:
759 :    
760 :     # Calculate equally weighted portfolio:
761 : wuertz 49 means = object$returns
762 :     covmat = object$cov
763 : wuertz 17 ew.weights = rep(1/length(object$returns), times = length(object$returns))
764 :     Rew = (ew.weights %*% means)[[1,1]]
765 :     Sew = (ew.weights %*% covmat %*% ew.weights)[[1,1]]
766 :     if (add) {
767 : wuertz 49 points(sqrt(Sew), Rew, col = "steelblue")
768 : wuertz 17 }
769 :    
770 :     # Return Value:
771 :     object$Rew = Rew
772 :     object$Sew = Sew
773 :     object
774 :     }
775 :    
776 :    
777 :     # ------------------------------------------------------------------------------
778 :     # Markowitz - Short Selling
779 :    
780 :    
781 :     .frontierShortSellingMarkowitz =
782 :     function(data, Rf = 0, length = 300,
783 :     r.range = NULL, s.range = NULL,
784 :     title = NULL, description = NULL, ...)
785 :     { # A function implemented by Diethelm Wuertz
786 :    
787 :     # Description:
788 :     # Calculates the efficient frontier from a matrix
789 :     # of either market or simulated assets given in matrix "x".
790 :     # Each time series represents a column in this matrix.
791 :    
792 :     # Arguments:
793 :     # data - either a rectangular time series object which can be
794 :     # transformed into a matrix or a list with two named entries,
795 :     # mu - the returns of the multivariate series, and
796 :     # Sigma - the Covariance matrix of the multivariate series
797 :     # length - Number of equidistant return points on the Efficient
798 :     # frontier
799 :     # r.range - Plot range of returns, if NULL, the range will be created
800 :     # automatically
801 :     # s.range - Plot range of standard deviations, if NULL, the range
802 :     # will be created automatically
803 :    
804 :     # FUNCTION:
805 :    
806 :     # Check Data:
807 :     if (is.list(data)) {
808 : wuertz 49 x.mat = NA
809 : wuertz 17 mu = data$mu
810 :     Sigma = data$Sigma
811 :     plottype = "mucov"
812 :     } else {
813 :     x.mat = as.matrix(x)
814 :     # Mean Vector and Covariance Matrix:
815 :     mu = apply(x.mat, MARGIN = 2, FUN = mean)
816 :     Sigma = apply(x.mat, MARGIN = 2, FUN = var)
817 :     plottype = "series"
818 :     }
819 :    
820 :     # Settings:
821 :     C0 = 1
822 :    
823 :     # Ranges for mean and Standard Deviation:
824 :     if (is.null(r.range)) r.range = range(mu)
825 :     if (is.null(s.range)) s.range = c(0, max(sqrt(diag(Sigma))))
826 :    
827 :     # Portfolio Settings:
828 :     one = rep(1, times = length(mu))
829 :     invSigma = solve(Sigma)
830 :     a = mu %*% invSigma %*% mu
831 :     b = mu %*% invSigma %*% one
832 :     c = one %*% invSigma %*% one
833 :     d = a*c - b^2
834 :    
835 :     # Minimum Variance Portfolio
836 :     mu.mv = (b/c)*C0
837 :     sigma.mv = C0/sqrt(c)
838 :    
839 :     # Plot:
840 :     mu.p = seq(r.range[1], r.range[2], length = length)
841 :     sigma.p = sqrt((c*mu.p^2 - 2*b*C0*mu.p + a*C0^2)/d)
842 :    
843 :     # Tangency Portfolio:
844 :     mu.tg = (a/b)*C0
845 :     sigma.tg = (sqrt(a)/b)*C0
846 :     weights.tg = C0 * as.vector(invSigma %*% mu ) / b
847 :    
848 :     # Market Portfolio:
849 :     mu.f = 0
850 :     mu.mk = 0
851 :     sigma.mk = 0
852 :    
853 :     # Equal Weights Portfolio:
854 :     ew.weights = rep(1/length(mu), times = length(mu))
855 :     mu.ew = (ew.weights %*% mu)[[1, 1]]
856 :     sigma.ew = (ew.weights %*% Sigma %*% ew.weights)[[1, 1]]
857 :    
858 :     # Result:
859 :     pfolio = list(
860 :     what = "Short Selling Markowitz",
861 :     plottype = plottype,
862 :     data = data,
863 :     pw = NA, pm = mu.p, ps = sigma.p,
864 :     returns = mu, cov = Sigma,
865 :     r.range = r.range, s.range = s.range,
866 :     # Tangency Portfolio:
867 :     Rf = Rf, Rm = mu.tg, Sm = sigma.tg,
868 :     # Equal Weights Portfolio:
869 :     Rew = mu.ew, Sew = sigma.ew,
870 :     t.weights = weights.tg,
871 :     diversification = NA)
872 :    
873 :     # Title:
874 :     if (is.null(title))
875 :     title = "Markowitz Portfolio Optimization"
876 :    
877 :     # Description:
878 :     if (is.null(description))
879 :     description = as.character(date())
880 :    
881 :     # Return Value:
882 :     new ("fPFOLIO",
883 :     call = as.call(match.call()),
884 :     method = "Analytical Solution",
885 :     model = "Short Selling Markowitz Portfolio",
886 :     data = list(x = data),
887 :     pfolio = pfolio,
888 :     title = as.character(title),
889 :     description = as.character(description)
890 :     )
891 :    
892 :     }
893 :    
894 :    
895 :     # ------------------------------------------------------------------------------
896 :    
897 :    
898 :     .frontierShortSellingMarkowitz.RUnit =
899 :     function()
900 :     { # A function implemented by Diethelm Wuertz
901 :    
902 :     # Description:
903 :     # RUnit Test
904 :    
905 :     # FUNCTION:
906 :    
907 :     # PortfolioNames:
908 :     Names = c("Elsevier", "Fortis", "Getronics", "Heineken", "Philips",
909 :     "RoyalDutch", "Unilever")
910 :    
911 :     # Mean:
912 :     mu = 1e-3 * c(0.266, 0.274, 0.162, 0.519, 0.394, 0.231, 0.277)
913 :     names(mu) = Names
914 :     print(mu)
915 :    
916 :     # Covariance Matrix:
917 :     Sigma = 1e-3 * matrix(c(
918 :     0.345, 0.150, 0.183, 0.088, 0.186, 0.090, 0.095,
919 :     0.150, 0.399, 0.204, 0.107, 0.236, 0.130, 0.127,
920 :     0.183, 0.204, 1.754, 0.075, 0.325, 0.110, 0.091,
921 :     0.088, 0.107, 0.075, 0.243, 0.096, 0.064, 0.086,
922 :     0.186, 0.236, 0.325, 0.096, 0.734, 0.147, 0.114,
923 :     0.090, 0.130, 0.110, 0.064, 0.147, 0.221, 0.093,
924 :     0.095, 0.127, 0.091, 0.086, 0.114, 0.093, 0.219),
925 :     byrow = TRUE, ncol = 7)
926 :     colnames(Sigma) = rownames(Sigma) = Names
927 :     print(Sigma)
928 :    
929 :     # Return Value:
930 :     .frontierShortSellingMarkowitz(data = list(mu = mu, Sigma = Sigma))
931 :     }
932 :    
933 :    
934 :     ################################################################################
935 :     # BUILTIN: tseries
936 :    
937 :    
938 :     .portfolio.optim =
939 : wuertz 49 function(pm, returns, covmat)
940 : wuertz 17 { # A Builtin function modified by Diethelm Wuertz
941 :    
942 :     # Description:
943 :     # Package: tseries
944 :     # Version: 0.9-21
945 :     # Date: 2004-04-23
946 :     # Title: Time series analysis and computational finance
947 :     # Author: Compiled by Adrian Trapletti <a.trapletti@bluewin.ch>
948 :     # Maintainer: Kurt Hornik <Kurt.Hornik@R-project.org>
949 :     # Description: Package for time series analysis and computational
950 :     # finance
951 :     # Depends: R (>= 1.9.0), quadprog
952 :     # License: GPL (see file COPYING)
953 :     # Packaged: Thu Apr 22 16:32:16 2004; hornik
954 :     # Built: R 1.9.0; i386-pc-mingw32; 2004-04-22 17:49:17; windows
955 :    
956 :     # FUNCTION:
957 :    
958 :     # Optimize:
959 : wuertz 49 k = dim(covmat)[2] # dim(x)[2]
960 : wuertz 17 dvec = rep(0, k)
961 :     a1 = rep(1, k)
962 : wuertz 49 a2 = returns # apply(x, 2, mean)
963 : wuertz 17 a3 = matrix(0, k, k)
964 :     diag(a3) = 1
965 :     b3 = rep(0, k)
966 :     Amat = t(rbind(a1, a2, a3))
967 :     b0 = c(1, pm, b3)
968 :     res = .solve.QP(covmat, dvec, Amat, bvec = b0, meq = 2)
969 : wuertz 49 weights = res$solution
970 :     ps = sqrt(weights %*% covmat %*% weights)
971 : wuertz 17
972 :     # Result:
973 : wuertz 49 ans = list(pw = weights, pm = pm, ps = ps)
974 : wuertz 17
975 :     # Return value:
976 :     ans
977 :     }
978 :    
979 :    
980 :     ################################################################################
981 :     # BUILTIN: quadprog
982 :    
983 :    
984 :     .solve.QP =
985 :     function(Dmat, dvec, Amat, bvec, meq)
986 :     { # A Builtin function modified by Diethelm Wuertz
987 :    
988 :     # Description:
989 :     # Package: quadprog
990 :     # Version: 1.4-7
991 :     # Date: 2004-01-31
992 :     # Title: Functions to solve Quadratic Programming Problems.
993 :     # Author: S original by Berwin A. Turlach <berwin.turlach@anu.edu.au>
994 :     # R port by Andreas Weingessel <Andreas.Weingessel@ci.tuwien.ac.at>
995 :     # Maintainer: Andreas Weingessel <Andreas.Weingessel@ci.tuwien.ac.at>
996 :     # Description: This package contains routines and documentation for
997 :     # solving quadratic programming problems.
998 :     # License: GPL version 2 or later
999 :     # Packaged: Sat Jan 31 13:32:53 2004; hornik
1000 :     # Built: R 1.9.0; i386-pc-mingw32; 2004-03-28 15:03:03
1001 :    
1002 :     # FUNCTION:
1003 :    
1004 :     # Solve QP:
1005 :     n = nrow(Dmat)
1006 :     q = ncol(Amat)
1007 :     iact = rep(0,q)
1008 :     nact = 0
1009 :     r = min(n,q)
1010 :     sol = rep(0,n)
1011 :     crval = 0
1012 :     work = rep(0,2*n+r*(r+5)/2+2*q+1)
1013 :     iter = rep(0,2)
1014 :     factorized = FALSE
1015 :     res1 = .Fortran("qpgen2",
1016 :     as.double(Dmat), dvec = as.double(dvec), as.integer(n),
1017 :     as.integer(n), sol = as.double(sol), crval = as.double(crval),
1018 :     as.double(Amat), as.double(bvec), as.integer(n), as.integer(q),
1019 :     as.integer(meq), iact = as.integer(iact), nact = as.integer(nact),
1020 :     iter = as.integer(iter), work = as.double(work),
1021 :     ierr = as.integer(factorized), PACKAGE = "fPortfolio")
1022 :     if (res1$ierr == 1) {
1023 : wuertz 49 warning("constraints are inconsistent, no solution!")
1024 : wuertz 17 } else if (res1$ierr == 2) {
1025 : wuertz 49 warning("matrix D in quadratic function is not positive definite!")
1026 : wuertz 17 }
1027 :    
1028 :     # Return Value:
1029 :     list(solution = res1$sol, value = res1$crval,
1030 :     unconstrainted.solution = res1$dvec,
1031 :     iterations = res1$iter, iact = res1$iact[1:res1$nact])
1032 :     }
1033 :    
1034 :    
1035 :     ################################################################################
1036 :    

R-Forge@R-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business University of Wisconsin - Madison Powered By FusionForge