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 49 - (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 :     # .frontier.fPFOLIO S3: Extract points on the efficient frontier
38 :     # 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 :     return(.frontierShortSelling(data = data, Rf = Rf, length = length,
120 :     r.range = r.range, s.range = s.range, title = title,
121 :     description = description, ...))
122 :     }
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 :     .frontier.fPFOLIO =
363 :     function(object)
364 :     { # A function implemented by Diethelm Wuertz
365 :    
366 :     # Description:
367 :     # S3 method for an object of class "fPFOLIO" to extract
368 :     # the points on the efficient frontier
369 :    
370 :     # FUNCTION:
371 :    
372 :     # Frontier:
373 :     pfolio = object@pfolio
374 :     mu = pfolio$pm
375 :     Sigma = pfolio$ps
376 :    
377 :     # Return Value:
378 :     list(Sigma = Sigma, mu = mu)
379 :     }
380 :    
381 :    
382 :     # ******************************************************************************
383 :    
384 :    
385 :     print.fPFOLIO =
386 :     function(x, ...)
387 :     { # A function implemented by Diethelm Wuertz
388 :    
389 :     # Description:
390 :     # S3 Print Method for an object of class "fPFOLIO"
391 :    
392 :     # Arguments:
393 :     # x - an object of class "fPFOLIO"
394 :    
395 :     # FUNCTION:
396 :    
397 :     # Extract Portfolio:
398 :     pfolio = x@pfolio
399 :    
400 :     # Title:
401 :     cat("\nTitle:\n ")
402 :     cat(as.character(x@title), "\n")
403 :    
404 :     # Call:
405 :     cat("\nCall:\n ")
406 :     print.default(x@call)
407 :    
408 :     # Portfolio Optimization:
409 :     if (pfolio$what == "portfolio") {
410 :     # Weights:
411 :     weights = round(pfolio$pw, digits = 2)
412 :     names(weights) = as.character(1:length(weights))
413 :     cat("\nPortfolio Weights:\n ")
414 :     print(weights[weights > 0])
415 :     # Check Sum of Weights:
416 :     cat("\nSum of Weights:\n ")
417 :     print(round(sum(pfolio$pw), digits = 2))
418 :     # Target Returns:
419 :     cat("\nTarget Return(s):\n ")
420 :     print(round(pfolio$pm, digits = 4))
421 :     # Target Risk:
422 :     cat("\nTarget Risk(s):\n ")
423 :     print(round(pfolio$ps, digits = 4))
424 :     }
425 :    
426 :     # Frontier:
427 :     if (pfolio$what == "frontier") {
428 :     # Efficient Frontier:
429 :     cat("\nEfficient Frontier - Returns:\n ")
430 :     r.digits = abs(round(log10(max(pfolio$pm)), 0)) + 4
431 :     print(round(pfolio$pm, digits = r.digits))
432 :     cat("\nEfficient Frontier - Standard Deviations:\n ")
433 :     s.digits = abs(round(log10(max(pfolio$ps)), 0)) + 4
434 :     print(round(pfolio$ps, digits = s.digits))
435 :     # Tangency Portfolio:
436 :     cat("\nTangency Portfolio:\n ")
437 :     cat("Risk free rate:", round(pfolio$Rf, digits = r.digits),
438 :     " Return:", round(pfolio$Rm, digits = r.digits),
439 :     " Risk:", round(pfolio$Sm, digits = s.digits), "\n")
440 :     # Equal Weights Portfolio:
441 :     cat("\nEqual Weights Portfolio:\n ")
442 :     cat("Risk:", round(pfolio$Rew, digits = r.digits),
443 :     " Return:", round(pfolio$Sew, digits = s.digits), "\n")
444 :     }
445 :    
446 :     # Description:
447 :     cat("\nDescription:\n ")
448 :     print(x@description)
449 :    
450 :     # Return Value:
451 :     invisible(x)
452 :     }
453 :    
454 :    
455 :     # ------------------------------------------------------------------------------
456 :    
457 :    
458 :     plot.fPFOLIO =
459 :     function(x, alpha = 0.05, mc = 500, which = "ask", ...)
460 :     { # A function implemented by Diethelm Wuertz
461 :    
462 :     # Description:
463 :     # S3 Print Method for an object of class "fPFOLIO"
464 :    
465 :     # Arguments:
466 :     # x - an object of class "fPFOLIO"
467 :    
468 :     # FUNCTION:
469 :    
470 :     # Check:
471 :     if (as.character(x@call[1]) == "portfolioMarkowitz")
472 :     stop("Objects returned by portfolioMarkowitz() cannot be plotted")
473 :    
474 :     # Settings:
475 :     .mcSteps <<- mc
476 :     .alphaLevel <<- alpha
477 :     if (which[1] == "ask") {
478 :     .interactive <<- TRUE
479 :     } else {
480 :     .interactive <<- FALSE
481 :     }
482 :    
483 :     # Plot Efficient Frontier:
484 :     plot.1 <<- function(x) {
485 :     pfolio = x@pfolio
486 :     # Extract:
487 :     dim = length(pfolio$returns)
488 :     returns = pfolio$returns
489 :     covar = pfolio$cov
490 :     r.range = pfolio$r.range
491 :     s.range = pfolio$s.range
492 :     diversification = pfolio$diversification
493 :     # 1. Plot Graph Frame:
494 :     plot(x = s.range, y = r.range, type = "n",
495 :     xlab = "Standard Deviation", ylab = "Return",
496 :     main = "Mean Variance Portfolio")
497 :     # 2. Plot Individual Assets:
498 :     points(sqrt(diag(covar)), returns, pch = 20, col = "red")
499 :     risk = sqrt(diag(covar)) + diff(s.range) / 50
500 :     text(risk, returns, as.character(1:dim), col = "red")
501 :     # 3. Plot Efficient Frontier:
502 :     points(pfolio$ps, pfolio$pm, col = "blue", cex = 0.1)
503 :     # 4. Add Tangency Portfolio:
504 :     # Find the Index (location of the Tangency Point):
505 :     delta.pm = diff(pfolio$pm)
506 :     delta.ps = diff(pfolio$ps)
507 :     slope1 = slope.keep = delta.pm/delta.ps
508 :     slope1 = slope1[slope.keep > 0]
509 :     slope2 = (pfolio$pm-pfolio$Rf) / pfolio$ps
510 :     diff.slope2 = diff(slope2) / 2
511 :     slope2 = slope2[1:length(diff.slope2)] + diff.slope2
512 :     slope2 = slope2[slope.keep > 0]
513 :     cross = abs(slope2 - slope1)
514 :     index = order(cross)[1]
515 :     pm.pos = pfolio$pm[slope.keep > 0]
516 :     ps.pos = pfolio$ps[slope.keep > 0]
517 :     # Here are the mean return and the standard deviation
518 :     mean.m = (pm.pos[index] + pm.pos[index+1]) / 2
519 :     sigma.m = (ps.pos[index] + ps.pos[index+1]) / 2
520 :     # ... together with the slope
521 :     slope = (mean.m - pfolio$Rf) / sigma.m
522 :     # Add the Capital Market Line and the Tangency Point to the plot
523 :     lines(c(0, max(pfolio$ps)), c(pfolio$Rf, slope* max(pfolio$ps) +
524 :     pfolio$Rf), lwd = 2 )
525 :     points(sigma.m, slope*sigma.m + pfolio$Rf, col = "green",
526 :     pch = 20, cex = 1.5)
527 :     s.delta = max(pfolio$ps) / 10
528 :     text(sigma.m - s.delta, slope*sigma.m + pfolio$Rf, "CML")
529 :     # 5. Add Equal Weights Portfolio:
530 :     # Calculate equally weighted portfolio:
531 :     ew.weights = rep(1/dim, times = dim)
532 :     Rew = (ew.weights %*% returns)
533 :     Sew = (ew.weights %*% covar %*% ew.weights)
534 :     points(sqrt(Sew), Rew, pch = 20, col = "steelblue4", cex = 1.5)
535 :     }
536 :     # ... add Monte Carlo Portfolios:
537 :     plot.2 <<- function(x) {
538 :     plot.1(x)
539 :     if (.interactive) {
540 :     mc = readline("Number of Monte Carlo Steps > ")
541 :     if (mc != "") .mcSteps <<- as.integer(mc)
542 :     cat("\nTotal Number of Monte Carlo Steps:", .mcSteps, "\n")
543 :     }
544 :     tmp = montecarloMarkowitz(x, mc = .mcSteps, doplot = FALSE, add = TRUE)
545 :    
546 :     }
547 :     # Bar Plot of Tangency Weights:
548 :     plot.3 <<- function(x) {
549 :     pfolio = x@pfolio
550 :     weights = pfolio$t.weights ## [pfolio$t.weights > 0]
551 :     ## pie(weights, col = rainbow(length(weights)), radius = 0.9,
552 :     ## main = "Tangency Portfolio Weights")
553 :     barplot(weights, cex.names = 0.7, col = "steelblue")
554 :     abline (h = 0)
555 :     }
556 :     # Plot Diversification:
557 :     plot.4 <<- function(x) {
558 :     pfolio = x@pfolio
559 :     if (is.na(pfolio$diversification)) {
560 :     cat ("\nNo time series data available.\n")
561 :     } else {
562 :     xx = pfolio$pm
563 :     yy = pfolio$diversification
564 :     data = data.frame(cbind(xx, yy))
565 :     plot(pfolio$pm, pfolio$diversification, type = "l", xlab = "Return",
566 :     ylab = "Number of Assets", main = "Diversification")
567 :     zz = loess(yy ~ xx, data)
568 :     lines(xx, zz$fit, col = "blue")
569 :     points(pfolio$Rm, length(weights), col = "green", cex = 1.5)
570 :     }
571 :     }
572 :     # 5. Plot Cumulated Return:
573 :     plot.5 <<- function(x) {
574 :     pfolio = x@pfolio
575 :     if (pfolio$plottype == "mucov") {
576 :     cat ("\nNo time series data available.\n")
577 :     } else {
578 :     plot(cumsum(as.matrix(pfolio$data) %*% pfolio$t.weights),
579 :     type = "l", col = "steelblue",
580 :     ylab = "Cumulated Return", main = "Portfolio Return")
581 :     grid()
582 :     }
583 :     }
584 :     # 6. Portfolio Risk Histogram:
585 :     plot.6 <<- function(x) {
586 :     pfolio = x@pfolio
587 :     if (pfolio$plottype == "mucov") {
588 :     cat ("\nNo time series data available.\n")
589 :     } else {
590 :     weights = pfolio$t.weights[pfolio$t.weights > 0]
591 :     pfolioHist(x@data, weights = weights, alpha = .alphaLevel)
592 :     }
593 :     }
594 :    
595 :     # Plot:
596 :     interactivePlot(
597 :     x = x,
598 :     choices = c(
599 :     "Efficient Frontier",
600 :     "... with Monte Carlo Portfolios",
601 :     "Bar Chart of Tangency Portfolio",
602 :     "Diversification Plot",
603 :     "Cumulated Return Tangency Portfolio",
604 :     "Risk of Tangency Portfolio"),
605 :     plotFUN = c(
606 :     "plot.1",
607 :     "plot.2",
608 :     "plot.3",
609 :     "plot.4",
610 :     "plot.5",
611 :     "plot.6"),
612 :     which = which)
613 :    
614 :    
615 :     # Return Value:
616 :     invisible(x)
617 :     }
618 :    
619 :    
620 :     # ------------------------------------------------------------------------------
621 :    
622 :    
623 :     summary.fPFOLIO =
624 :     function(object, ...)
625 :     { # A function implemented by Diethelm Wuertz
626 :    
627 :     # Description:
628 :     # Summary Method for an object of class "fPFOLIO"
629 :    
630 :     # FUNCTION:
631 :    
632 :     # Extract Portfolio:
633 :     pfolio = object@pfolio
634 :    
635 :     # Print Result:
636 :     print(object)
637 :    
638 :     # Plot:
639 :     plot(object, ...)
640 :    
641 :     # Return Value:
642 :     invisible(object)
643 :     }
644 :    
645 :    
646 :     ################################################################################
647 :    
648 :    
649 :     .tangencyMarkowitz =
650 :     function(object, Rf = 0, add = TRUE)
651 :     { # A function implemented by Diethelm Wuertz
652 :    
653 :     # Description:
654 :     # This function calculates the properties of the tangency portfolio,
655 :     # and adds the "capital market line" and the tangency portfolio, Rm
656 :     # and Sm, to the efficient frontier.
657 :    
658 :     # NOTE:
659 :     # Not for all risk free rates there exists a tangency portfolio!!!
660 :    
661 :     # Arguments:
662 :     # object - An object of class portfolio
663 :    
664 :     # Value:
665 :     # An updated object of class portfolio with the following
666 :     # components added/updated:
667 :     # Rf - risk free rate of return
668 :     # Rm - rate of return of the tangency portfolio
669 :     # Sm - standardard deviation of the tangency portfolio
670 :     # t.weights - weights of the tangency portfolio
671 :    
672 :     # Note:
673 :     # I have not yet checked if this works for all values of Rf !!
674 :     # I hope it works properly even if there exists no tangency portfolio.
675 :    
676 :     # FUNCTION:
677 :    
678 :     # Find the Index (location of the Tangency Point:
679 :     delta.pm = diff(object$pm)
680 :     delta.ps = diff(object$ps)
681 :     slope1 = slope.keep = delta.pm/delta.ps
682 :     slope1 = slope1[slope.keep > 0]
683 :     slope2 = (object$pm-Rf) / object$ps
684 :     diff.slope2 = diff(slope2) / 2
685 :     slope2 = slope2[1:length(diff.slope2)] + diff.slope2
686 :     slope2 = slope2[slope.keep > 0]
687 :     cross = abs(slope2 - slope1)
688 :     index = order(cross)[1]
689 :     pm.pos = object$pm[slope.keep > 0]
690 :     ps.pos = object$ps[slope.keep > 0]
691 :    
692 :     # Here are the mean return and the standard deviation
693 :     mean.m = (pm.pos[index] + pm.pos[index+1]) / 2
694 :     sigma.m = (ps.pos[index] + ps.pos[index+1]) / 2
695 :     # ... together with the slope
696 :     slope = (mean.m - Rf) / sigma.m
697 :     x1 = 0
698 :     y1 = Rf
699 :     x2 = max(object$ps)
700 :     y2 = slope*x2 + Rf
701 :     xm = sigma.m
702 :     ym = slope*sigma.m + Rf
703 :    
704 :     # Add the Capital Market Line and the Tangency Point to the plot
705 :     if (add) {
706 :     lines(c(x1, x2), c(y1, y2), lwd = 2 )
707 :     points(xm, ym, col = "green", cex = 1.5)
708 :     s.delta = max(object$ps) / 10
709 :     text(xm - s.delta, ym, "CML")}
710 :    
711 :     # Calculate Tangency Portfolio, if it exists!
712 :     # The test is true, when the slope could be evaluated.
713 :     if (!is.na(slope)) {
714 : wuertz 49 tangency = .portfolio.optim(pm = ym, returns = object$returns,
715 : wuertz 17 covmat = object$cov)
716 :     # Portfolio Weights in Percent:
717 :     pw = 100 * tangency$pw
718 :     # Remove those smaller 1 %% Promille!!
719 :     pw = round(pw*(sign(pw - 0.1) + 1)/2, digits = 2)
720 :     object$t.weights = pw/100
721 :     }
722 :    
723 :     # Add to Output:
724 :     object$Rf = Rf
725 :     object$Rm = ym
726 :     object$Sm = xm
727 :    
728 :     # Return Value:
729 :     object
730 :     }
731 :    
732 :    
733 :     # ------------------------------------------------------------------------------
734 :    
735 :    
736 :     .equalweightsMarkowitz =
737 :     function(object, add = TRUE)
738 :     { # A function implemented by Diethelm Wuertz
739 :    
740 :     # Description:
741 :     # This function adds the equal weighted Portfolio to
742 :     # efficient frontier plot.
743 :    
744 :     # FUNCTION:
745 :    
746 :     # Calculate equally weighted portfolio:
747 : wuertz 49 means = object$returns
748 :     covmat = object$cov
749 : wuertz 17 ew.weights = rep(1/length(object$returns), times = length(object$returns))
750 :     Rew = (ew.weights %*% means)[[1,1]]
751 :     Sew = (ew.weights %*% covmat %*% ew.weights)[[1,1]]
752 :     if (add) {
753 : wuertz 49 points(sqrt(Sew), Rew, col = "steelblue")
754 : wuertz 17 }
755 :    
756 :     # Return Value:
757 :     object$Rew = Rew
758 :     object$Sew = Sew
759 :     object
760 :     }
761 :    
762 :    
763 :     # ------------------------------------------------------------------------------
764 :     # Markowitz - Short Selling
765 :    
766 :    
767 :     .frontierShortSellingMarkowitz =
768 :     function(data, Rf = 0, length = 300,
769 :     r.range = NULL, s.range = NULL,
770 :     title = NULL, description = NULL, ...)
771 :     { # A function implemented by Diethelm Wuertz
772 :    
773 :     # Description:
774 :     # Calculates the efficient frontier from a matrix
775 :     # of either market or simulated assets given in matrix "x".
776 :     # Each time series represents a column in this matrix.
777 :    
778 :     # Arguments:
779 :     # data - either a rectangular time series object which can be
780 :     # transformed into a matrix or a list with two named entries,
781 :     # mu - the returns of the multivariate series, and
782 :     # Sigma - the Covariance matrix of the multivariate series
783 :     # length - Number of equidistant return points on the Efficient
784 :     # frontier
785 :     # r.range - Plot range of returns, if NULL, the range will be created
786 :     # automatically
787 :     # s.range - Plot range of standard deviations, if NULL, the range
788 :     # will be created automatically
789 :    
790 :     # FUNCTION:
791 :    
792 :     # Check Data:
793 :     if (is.list(data)) {
794 : wuertz 49 x.mat = NA
795 : wuertz 17 mu = data$mu
796 :     Sigma = data$Sigma
797 :     plottype = "mucov"
798 :     } else {
799 :     x.mat = as.matrix(x)
800 :     # Mean Vector and Covariance Matrix:
801 :     mu = apply(x.mat, MARGIN = 2, FUN = mean)
802 :     Sigma = apply(x.mat, MARGIN = 2, FUN = var)
803 :     plottype = "series"
804 :     }
805 :    
806 :     # Settings:
807 :     C0 = 1
808 :    
809 :     # Ranges for mean and Standard Deviation:
810 :     if (is.null(r.range)) r.range = range(mu)
811 :     if (is.null(s.range)) s.range = c(0, max(sqrt(diag(Sigma))))
812 :    
813 :     # Portfolio Settings:
814 :     one = rep(1, times = length(mu))
815 :     invSigma = solve(Sigma)
816 :     a = mu %*% invSigma %*% mu
817 :     b = mu %*% invSigma %*% one
818 :     c = one %*% invSigma %*% one
819 :     d = a*c - b^2
820 :    
821 :     # Minimum Variance Portfolio
822 :     mu.mv = (b/c)*C0
823 :     sigma.mv = C0/sqrt(c)
824 :    
825 :     # Plot:
826 :     mu.p = seq(r.range[1], r.range[2], length = length)
827 :     sigma.p = sqrt((c*mu.p^2 - 2*b*C0*mu.p + a*C0^2)/d)
828 :    
829 :     # Tangency Portfolio:
830 :     mu.tg = (a/b)*C0
831 :     sigma.tg = (sqrt(a)/b)*C0
832 :     weights.tg = C0 * as.vector(invSigma %*% mu ) / b
833 :    
834 :     # Market Portfolio:
835 :     mu.f = 0
836 :     mu.mk = 0
837 :     sigma.mk = 0
838 :    
839 :     # Equal Weights Portfolio:
840 :     ew.weights = rep(1/length(mu), times = length(mu))
841 :     mu.ew = (ew.weights %*% mu)[[1, 1]]
842 :     sigma.ew = (ew.weights %*% Sigma %*% ew.weights)[[1, 1]]
843 :    
844 :     # Result:
845 :     pfolio = list(
846 :     what = "Short Selling Markowitz",
847 :     plottype = plottype,
848 :     data = data,
849 :     pw = NA, pm = mu.p, ps = sigma.p,
850 :     returns = mu, cov = Sigma,
851 :     r.range = r.range, s.range = s.range,
852 :     # Tangency Portfolio:
853 :     Rf = Rf, Rm = mu.tg, Sm = sigma.tg,
854 :     # Equal Weights Portfolio:
855 :     Rew = mu.ew, Sew = sigma.ew,
856 :     t.weights = weights.tg,
857 :     diversification = NA)
858 :    
859 :     # Title:
860 :     if (is.null(title))
861 :     title = "Markowitz Portfolio Optimization"
862 :    
863 :     # Description:
864 :     if (is.null(description))
865 :     description = as.character(date())
866 :    
867 :     # Return Value:
868 :     new ("fPFOLIO",
869 :     call = as.call(match.call()),
870 :     method = "Analytical Solution",
871 :     model = "Short Selling Markowitz Portfolio",
872 :     data = list(x = data),
873 :     pfolio = pfolio,
874 :     title = as.character(title),
875 :     description = as.character(description)
876 :     )
877 :    
878 :     }
879 :    
880 :    
881 :     # ------------------------------------------------------------------------------
882 :    
883 :    
884 :     .frontierShortSellingMarkowitz.RUnit =
885 :     function()
886 :     { # A function implemented by Diethelm Wuertz
887 :    
888 :     # Description:
889 :     # RUnit Test
890 :    
891 :     # FUNCTION:
892 :    
893 :     # PortfolioNames:
894 :     Names = c("Elsevier", "Fortis", "Getronics", "Heineken", "Philips",
895 :     "RoyalDutch", "Unilever")
896 :    
897 :     # Mean:
898 :     mu = 1e-3 * c(0.266, 0.274, 0.162, 0.519, 0.394, 0.231, 0.277)
899 :     names(mu) = Names
900 :     print(mu)
901 :    
902 :     # Covariance Matrix:
903 :     Sigma = 1e-3 * matrix(c(
904 :     0.345, 0.150, 0.183, 0.088, 0.186, 0.090, 0.095,
905 :     0.150, 0.399, 0.204, 0.107, 0.236, 0.130, 0.127,
906 :     0.183, 0.204, 1.754, 0.075, 0.325, 0.110, 0.091,
907 :     0.088, 0.107, 0.075, 0.243, 0.096, 0.064, 0.086,
908 :     0.186, 0.236, 0.325, 0.096, 0.734, 0.147, 0.114,
909 :     0.090, 0.130, 0.110, 0.064, 0.147, 0.221, 0.093,
910 :     0.095, 0.127, 0.091, 0.086, 0.114, 0.093, 0.219),
911 :     byrow = TRUE, ncol = 7)
912 :     colnames(Sigma) = rownames(Sigma) = Names
913 :     print(Sigma)
914 :    
915 :     # Return Value:
916 :     .frontierShortSellingMarkowitz(data = list(mu = mu, Sigma = Sigma))
917 :     }
918 :    
919 :    
920 :     ################################################################################
921 :     # BUILTIN: tseries
922 :    
923 :    
924 :     .portfolio.optim =
925 : wuertz 49 function(pm, returns, covmat)
926 : wuertz 17 { # A Builtin function modified by Diethelm Wuertz
927 :    
928 :     # Description:
929 :     # Package: tseries
930 :     # Version: 0.9-21
931 :     # Date: 2004-04-23
932 :     # Title: Time series analysis and computational finance
933 :     # Author: Compiled by Adrian Trapletti <a.trapletti@bluewin.ch>
934 :     # Maintainer: Kurt Hornik <Kurt.Hornik@R-project.org>
935 :     # Description: Package for time series analysis and computational
936 :     # finance
937 :     # Depends: R (>= 1.9.0), quadprog
938 :     # License: GPL (see file COPYING)
939 :     # Packaged: Thu Apr 22 16:32:16 2004; hornik
940 :     # Built: R 1.9.0; i386-pc-mingw32; 2004-04-22 17:49:17; windows
941 :    
942 :     # FUNCTION:
943 :    
944 :     # Optimize:
945 : wuertz 49 k = dim(covmat)[2] # dim(x)[2]
946 : wuertz 17 dvec = rep(0, k)
947 :     a1 = rep(1, k)
948 : wuertz 49 a2 = returns # apply(x, 2, mean)
949 : wuertz 17 a3 = matrix(0, k, k)
950 :     diag(a3) = 1
951 :     b3 = rep(0, k)
952 :     Amat = t(rbind(a1, a2, a3))
953 :     b0 = c(1, pm, b3)
954 :     res = .solve.QP(covmat, dvec, Amat, bvec = b0, meq = 2)
955 : wuertz 49 weights = res$solution
956 :     ps = sqrt(weights %*% covmat %*% weights)
957 : wuertz 17
958 :     # Result:
959 : wuertz 49 ans = list(pw = weights, pm = pm, ps = ps)
960 : wuertz 17
961 :     # Return value:
962 :     ans
963 :     }
964 :    
965 :    
966 :     ################################################################################
967 :     # BUILTIN: quadprog
968 :    
969 :    
970 :     .solve.QP =
971 :     function(Dmat, dvec, Amat, bvec, meq)
972 :     { # A Builtin function modified by Diethelm Wuertz
973 :    
974 :     # Description:
975 :     # Package: quadprog
976 :     # Version: 1.4-7
977 :     # Date: 2004-01-31
978 :     # Title: Functions to solve Quadratic Programming Problems.
979 :     # Author: S original by Berwin A. Turlach <berwin.turlach@anu.edu.au>
980 :     # R port by Andreas Weingessel <Andreas.Weingessel@ci.tuwien.ac.at>
981 :     # Maintainer: Andreas Weingessel <Andreas.Weingessel@ci.tuwien.ac.at>
982 :     # Description: This package contains routines and documentation for
983 :     # solving quadratic programming problems.
984 :     # License: GPL version 2 or later
985 :     # Packaged: Sat Jan 31 13:32:53 2004; hornik
986 :     # Built: R 1.9.0; i386-pc-mingw32; 2004-03-28 15:03:03
987 :    
988 :     # FUNCTION:
989 :    
990 :     # Solve QP:
991 :     n = nrow(Dmat)
992 :     q = ncol(Amat)
993 :     iact = rep(0,q)
994 :     nact = 0
995 :     r = min(n,q)
996 :     sol = rep(0,n)
997 :     crval = 0
998 :     work = rep(0,2*n+r*(r+5)/2+2*q+1)
999 :     iter = rep(0,2)
1000 :     factorized = FALSE
1001 :     res1 = .Fortran("qpgen2",
1002 :     as.double(Dmat), dvec = as.double(dvec), as.integer(n),
1003 :     as.integer(n), sol = as.double(sol), crval = as.double(crval),
1004 :     as.double(Amat), as.double(bvec), as.integer(n), as.integer(q),
1005 :     as.integer(meq), iact = as.integer(iact), nact = as.integer(nact),
1006 :     iter = as.integer(iter), work = as.double(work),
1007 :     ierr = as.integer(factorized), PACKAGE = "fPortfolio")
1008 :     if (res1$ierr == 1) {
1009 : wuertz 49 warning("constraints are inconsistent, no solution!")
1010 : wuertz 17 } else if (res1$ierr == 2) {
1011 : wuertz 49 warning("matrix D in quadratic function is not positive definite!")
1012 : wuertz 17 }
1013 :    
1014 :     # Return Value:
1015 :     list(solution = res1$sol, value = res1$crval,
1016 :     unconstrainted.solution = res1$dvec,
1017 :     iterations = res1$iter, iact = res1$iact[1:res1$nact])
1018 :     }
1019 :    
1020 :    
1021 :     ################################################################################
1022 :    

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