diff -Nru fassets-2100.78/ChangeLog fassets-2110.79/ChangeLog --- fassets-2100.78/ChangeLog 2009-09-28 12:27:20.000000000 +0000 +++ fassets-2110.79/ChangeLog 2012-09-24 07:39:47.000000000 +0000 @@ -1,78 +1,189 @@ -2009-05-06 wuertz +2012-09-24 chalabi - * shrink from tawny added als builtin + * DESCRIPTION: Updated maintainer field. + * NAMESPACE, R/zzz.Deprecated.R: Removed external C call. -pkg/fAssets/R/builtin-shrinkTawny.R +2011-09-23 mmaechler -2009-04-29 wuertz + * DESCRIPTION: remove deprecated "LazyLoad" entry - * comment added +2010-10-26 chalabi -pkg/fAssets/R/assetsMeanCov.R + * NAMESPACE: updated NAMESPACE -2009-04-19 chalabi +2010-07-23 chalabi - * added explicit version number in Depends field for key packages + * inst/DocCopying.pdf: removed DocCopying.pdf license is already + specified in DESCRIPTION file -pkg/fAssets/DESCRIPTION +2010-04-14 chalabi -2009-04-08 ellis + * NAMESPACE: updated NAMESPACE - * added function to compute color space for correlation matrix plot +2010-04-11 wuertz -pkg/fAssets/R/plot-pairs.R + * R/assets-arrange.R, R/assets-fit.R, R/assets-lpm.R, + R/assets-mcr.R, R/assets-meancov.R, R/assets-outliers.R, + R/assets-portfolio.R, R/assets-select.R, R/assets-simulate.R, + R/assets-test.R, R/assetsArrange.R, R/assetsFit.R, R/assetsLPM.R, + R/assetsMCR.R, R/assetsMeanCov.R, R/assetsOutliers.R, + R/assetsPfolio.R, R/assetsSelect.R, R/assetsSim.R, + R/assetsTest.R: some files renamed for consistency -2009-03-13 wuertz +2010-02-21 dscott - * small fix done + * ChangeLog, R/zzz.R, src/Makevars: Minor changes so passes check + * ChangeLog, R/plot-pairs.R, R/zzz.R, src/Makevars: minor changes + to plot-pairs.R and fixed dll with alteration to zzz.R and + src/Makevars -pkg/fAssets/R/assetsTest.R +2009-11-22 wuertz -2009-02-08 wuertz + * R/stats-distance.R: mutinfo function modified + * R/zzz.R, src/ecodist.c: C code modified to work all C programs + together + * R/stats-distance.R, src/ecodist.c: C Code added ecodist.c + * NAMESPACE, R/builtin-ecodist.R, R/stats-distance.R: namespace + updated + * R/builtin-ecodist.R: function name modified + * R/stats-distance.R: code modified + * R/builtin-ecodist.R, R/stats-distance.R: distance measures added + (undocumented) - * package reorganized, script files +2009-10-26 wuertz -pkg/fAssets/R/builtin-corrgram.R -pkg/fAssets/R/buitin-corrgram.R + * man/assetsMCR.Rd: man page examples corrected + * NAMESPACE, R/assetsMCR.R, R/assetsPfolio.R, R/zzz.Deprecated.R, + man/assetsMCR.Rd: assetsMCR.R and .Rd script added for marginal + contribution to covariance risk -2009-02-08 wuertz +2009-09-28 chalabi - * script files reorganized + * DESCRIPTION: updated version number + * ChangeLog, DESCRIPTION: updated DESCR and ChangeLog + * NAMESPACE: new NAMESPACE structure which should ease maintenance + of packages. + * DESCRIPTION, NAMESPACE: Merge branch 'devel-timeSeries' + + Conflicts: + pkg/timeSeries/R/base-Extract.R + pkg/timeSeries/R/timeSeries.R -pkg/fAssets/R/assetsArrange.R -pkg/fAssets/R/assetsFit.R -pkg/fAssets/R/assetsLPM.R -pkg/fAssets/R/assetsMeanCov.R -pkg/fAssets/R/assetsOutliers.R -pkg/fAssets/R/assetsPfolio.R -pkg/fAssets/R/assetsSelect.R -pkg/fAssets/R/assetsSim.R -pkg/fAssets/R/assetsStats.R -pkg/fAssets/R/assetsTest.R -pkg/fAssets/R/panel-diagonal.R -pkg/fAssets/R/plot-panels.R -pkg/fAssets/R/plot-qqplot.R -pkg/fAssets/R/zzz.Deprecated.R +2009-05-06 wuertz -2009-02-08 wuertz + * R/builtin-shrinkTawny.R: shrink from tawny added als builtin + * R/assetsMeanCov.R: Mean Cov functionality extended - * script files freorganisation +2009-05-01 wuertz -pkg/fAssets/R/assetsFit.R -pkg/fAssets/R/assetsPfolio.R -pkg/fAssets/R/assetsRisk.R -pkg/fAssets/R/class-fASSETS.R -pkg/fAssets/R/plot-binning.R -pkg/fAssets/R/plot-boxplot.R -pkg/fAssets/R/plot-ellipses.R -pkg/fAssets/R/plot-hist.R -pkg/fAssets/R/plot-mst.R -pkg/fAssets/R/plot-pairs.R -pkg/fAssets/R/plot-panels.R -pkg/fAssets/R/plotPanels.R + * R/assetsMeanCov.R, R/builtin-arwMvoutlier.R: hidden robut + covariance stimator .arwMeanCov added -2009-01-02 wuertz +2009-04-29 wuertz - * stars plot, plot argument corrected + * R/assetsMeanCov.R: comment added + +2009-04-28 wuertz + + * R/builtin-mst.R, R/builtin-robust.R, man/assetsMeanCov.Rd: new + robust cov estimators added + * NAMESPACE: namespace new functions added + * R/builtin-DEoptim.R, R/builtin-corpcor.R, R/builtin-corrgram.R, + R/builtin-covRobust.R, R/builtin-donostahRobust.R, + R/builtin-energy.R, R/builtin-mstApe.R, R/builtin-rmtTawny.R, + R/builtin-shrinkTawny.R: more information added to builtin + function, builtins added for shrink and rmt from tawny and bayes + stein from alexios + * R/assetsMeanCov.R: bayes stein, ledoit wolf, and rmt covariance + estimator added + +2009-04-19 chalabi + + * DESCRIPTION: added explicit version number in Depends field for + key packages + +2009-04-08 ellis + + * R/plot-pairs.R: added function to compute color space for + correlation matrix plot + +2009-04-02 chalabi + + * DESCRIPTION: more explicit depends and suggests field in DESC + file. + * DESCRIPTION: updated DESC file + +2009-03-13 wuertz + + * R/assetsTest.R: small fix done + * R/assetsTest.R: fixed + +2009-02-09 wuertz + + * NAMESPACE, R/assetsArrange.R, R/assetsMeanCov.R, R/assetsStats.R, + R/plot-risk.R, R/plot-stars.R, R/zzz.Deprecated.R, R/zzz.R, + inst/unitTests/runit.AssetsMeanCov.R, man/00fAssets-package.Rd, + man/VaRModelling.Rd, man/assetsArrange.Rd, man/assetsFit.Rd, + man/assetsLPM.Rd, man/assetsMeanCov.Rd, man/assetsOutliers.Rd, + man/assetsPfolio.Rd, man/assetsSelect.Rd, man/assetsSim.Rd, + man/assetsStats.Rd, man/assetsTest.Rd, man/boxPlot.Rd, + man/class-fASSETS.Rd, man/covEllipsesPlot.Rd, man/pairsPlot.Rd, + man/plot-binning.Rd, man/plot-boxplot.Rd, man/plot-ellipses.Rd, + man/plot-hist.Rd, man/plot-mst.Rd, man/plot-pairs.Rd, + man/plot-qqplot.Rd, man/plot-risk.Rd, man/plot-series.Rd, + man/plot-similarity.Rd, man/plot-stars.Rd, man/seriesPlot.Rd, + man/similarityPlot.Rd, man/starsPlot.Rd: help pages and + documentation essentiall improved, all functions, arguments and + retur5ned values should now be documented + +2009-02-08 wuertz + + * R/builtin-corrgram.R, R/buitin-corrgram.R: package reorganized, + script files + * R/assetsArrange.R, R/assetsFit.R, R/assetsLPM.R, + R/assetsMeanCov.R, R/assetsOutliers.R, R/assetsPfolio.R, + R/assetsSelect.R, R/assetsSim.R, R/assetsStats.R, R/assetsTest.R, + R/panel-diagonal.R, R/plot-panels.R, R/plot-qqplot.R, + R/zzz.Deprecated.R: script files reorganized + * R/panel-diagonal.R, R/plot-panels.R, R/plotPanels.R: + reorganization of files + * R/assetsFit.R, R/assetsPfolio.R, R/assetsRisk.R, + R/class-fASSETS.R, R/plot-binning.R, R/plot-boxplot.R, + R/plot-ellipses.R, R/plot-hist.R, R/plot-mst.R, R/plot-pairs.R, + R/plot-panels.R, R/plotPanels.R: script files freorganisation + * R/VaRModelling.R, R/assetsRisk.R, R/exampleCovData.R, + R/fixBinHistogram.R, R/plot-binning.R, R/plot-correlation.R, + R/plot-covEllipses.R, R/plot-ellipses.R, R/plot-histPairs.R, + R/plot-minSpanTree.R, R/plot-mst.R, R/plot-pairsPanels.R, + R/plot-panels.R: files renamed + * R/assetsArrange.R, R/assetsFit.R, R/assetsMeanCov.R, + R/assetsOutliers.R, R/assetsSelect.R, R/assetsStats.R, + R/assetsTest.R, R/builtin-robust.R, R/exampleCovData.R, + R/fixBinHistogram.R, R/plot-correlation.R, R/plot-histPairs.R, + R/plot-stars.R, R/zzz.Deprecated.R: reorginization of files + +2009-01-27 wuertz + + * R/plot-pairs.R: warnings hidden for pairs() if tick = 0 + +2009-01-16 chalabi + + * man/assetsLPM.Rd, man/assetsMeanCov.Rd, man/seriesPlot.Rd: fixed + warning with new Rd parser + +2009-01-04 wuertz + + * NAMESPACE, R/assetsMeanCov.R, R/assetsOutliers.R, + R/outlierDetection.R, man/assetsMeanCov.Rd, + man/assetsOutliers.Rd: internal function .assetsOutlierDetection + moved to assetsOutlier, documented and added to NAMESPACE + +2009-01-02 wuertz + + * R/plot-boxplot.R: default abline removed from box plot + * R/plot-stars.R: stars plot, plot argument corrected + +2008-12-31 wuertz + + * R/assetsSelect.R, man/assetsSelect.Rd: small modifications -pkg/fAssets/R/plot-stars.R diff -Nru fassets-2100.78/DESCRIPTION fassets-2110.79/DESCRIPTION --- fassets-2100.78/DESCRIPTION 2009-09-28 14:09:45.000000000 +0000 +++ fassets-2110.79/DESCRIPTION 2012-09-24 08:31:47.000000000 +0000 @@ -1,22 +1,21 @@ Package: fAssets -Version: 2100.78 -Revision: 4285 -Date: 2009-09-28 +Version: 2110.79 +Revision: 5367 +Date: 2012-09-24 Title: Rmetrics - Assets Selection and Modelling Author: Diethelm Wuertz and many others, see the SOURCE file Depends: R (>= 2.6.0), methods, sn, MASS, robustbase, timeDate, timeSeries, fBasics, fCopulae (>= 2100.77) Suggests: RUnit -Maintainer: Rmetrics Core Team +Maintainer: Yohan Chalabi Description: Environment for teaching "Financial Engineering and Computational Finance" NOTE: SEVERAL PARTS ARE STILL PRELIMINARY AND MAY BE CHANGED IN THE FUTURE. THIS TYPICALLY INCLUDES FUNCTION AND ARGUMENT NAMES, AS WELL AS DEFAULTS FOR ARGUMENTS AND RETURN VALUES. -LazyLoad: yes LazyData: yes License: GPL (>= 2) URL: http://www.rmetrics.org -Packaged: 2009-09-28 12:27:29 UTC; yankee +Packaged: 2012-09-24 07:39:48 UTC; yankee Repository: CRAN -Date/Publication: 2009-09-28 14:09:45 +Date/Publication: 2012-09-24 08:31:47 diff -Nru fassets-2100.78/MD5 fassets-2110.79/MD5 --- fassets-2100.78/MD5 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/MD5 2012-09-24 08:31:47.000000000 +0000 @@ -0,0 +1,76 @@ +6f76cf88fea9fabd3ca573729596fe98 *ChangeLog +b532c8ac33be92fe88e03dacf20db7e5 *DESCRIPTION +b033a7b3f6deac4b50f65333297566bc *NAMESPACE +21ca3984f8736b38581ea8e910ade5fb *R/assets-arrange.R +9dbed515b79c57ef6591e2c4ae112d6b *R/assets-fit.R +0cc744f0b627b20f7c4b8ef3893323b2 *R/assets-lpm.R +1cc90c11d4c31b171d91096e80fdbd9f *R/assets-mcr.R +fc137624c0534df616f1ab7f543fadb5 *R/assets-meancov.R +86a4a54059922fab478f002fc9866079 *R/assets-outliers.R +56df94c1a06b03961a37a7ea7e6e073c *R/assets-portfolio.R +7ecc26d8d6a160ad62d6797a95072b4c *R/assets-select.R +63c107670a288d666aac1e9832822ee4 *R/assets-simulate.R +0c35719e3673be4a2728754338bd6fbf *R/assets-test.R +31963962bea8a207ba501f123c04a0e2 *R/builtin-DEoptim.R +f573257a8fc7e6eb3ebe8b64eacc0ce2 *R/builtin-arwMvoutlier.R +5246456a49bc0834d905113ada7d207b *R/builtin-corpcor.R +8b9c08b997e7e0d830ab7296cbc8eea4 *R/builtin-corrgram.R +ac6299cdfe8467717028c284c69da265 *R/builtin-covRobust.R +9889bc2ba9be634a2b9f29f3f9a682fb *R/builtin-donostahRobust.R +60b5bbdb35369e05dd055402faafe397 *R/builtin-ecodist.R +eac428f89af108d95b02d670a729c5fa *R/builtin-energy.R +2bdfd75a42e5b859192efe652be34239 *R/builtin-mstApe.R +5b4262591b2e4419a100a0a8ba39a029 *R/builtin-rmtTawny.R +a0406ec5b5409bb44dd52a1b84c9db58 *R/builtin-shrinkTawny.R +12c33103982050308a93a9f459b2e71e *R/class-fASSETS.R +5eb9c6e068549c54803862b54a928e05 *R/plot-binning.R +bab842e1a7eef3c4327a83c91743b6fe *R/plot-boxplot.R +babbd395b31ed9ba83c3acf6b1336614 *R/plot-ellipses.R +94685adbd1dcb83f4483b64bc3831a38 *R/plot-hist.R +19003188473f404343e50495fb1824f9 *R/plot-mst.R +e78a8f8420c1b1d1292155fd6904fd83 *R/plot-pairs.R +4064e22c8c4cd03c96143a09ebf1d863 *R/plot-panels.R +a5934ed8b6321b4c5a83bec55322563e *R/plot-qqplot.R +b0970b505f1e0b15e1f6208cd9ff06fc *R/plot-risk.R +7e39197b96d538c1c44fbdcab39e7773 *R/plot-series.R +1e181f17f9bf22e59402bbd3153632eb *R/plot-similarity.R +e7af5f1290411e8c3a6f30e66db0d0af *R/plot-stars.R +85854f76a2caaec474a7dcccae508ae3 *R/stats-distance.R +489dd5049c648353724393e7dd6f2ab3 *R/zzz.Deprecated.R +a7dc77317fdb084a267697662890a6ac *R/zzz.R +6042b9c5e5bec3ecc1b6959cd2858b64 *inst/COPYRIGHT.html +7d215d06e1f0f0794a63a29b7a8b7935 *inst/unitTests/Makefile +d11a67f141f4eef8cb4d3f24d66bd5d6 *inst/unitTests/runTests.R +2bd25e3404bc4640cbd56d5fade0899c *inst/unitTests/runit.AssetsFit.R +d8e6cb81482240804a1993f789894d65 *inst/unitTests/runit.AssetsMeanCov.R +86d6d79a41057e41f34bcc184e2476de *inst/unitTests/runit.AssetsPlots.R +f4f0d316b9b71062619c5e089d41b798 *inst/unitTests/runit.AssetsSelect.R +1fe17670137e36b52e7a9e9caf46a021 *inst/unitTests/runit.AssetsTests.R +57afd3f50fe97e70d3398074a2cf408c *inst/unitTests/runit.LowerPartialMoments.R +fd9a8572ff1bf7c579828811f391572d *man/00fAssets-package.Rd +3452cc11df57dae3bc32fc405c496b2c *man/assetsArrange.Rd +4a32199421f9dec71a3b41370724fc6f *man/assetsFit.Rd +62f9da699232821ded513b1270d63633 *man/assetsLPM.Rd +eda0293c9d377ce56f0c1081d75fa15f *man/assetsMCR.Rd +ed16cf26a8275a3f0d223716444959c4 *man/assetsMeanCov.Rd +3e5a8fcb45666327ff0d8cf33800b822 *man/assetsOutliers.Rd +929fbb072be2ef58526467eb4d4f9b90 *man/assetsPfolio.Rd +4e35a41e88825a646c8899f6d4286c54 *man/assetsSelect.Rd +52450aa784ba295d9c138b716813fae8 *man/assetsSim.Rd +0cad2b9cf5f5c210820732cb23cf20d4 *man/assetsTest.Rd +95b044331ae5e1c566a2126df92d6dcc *man/class-fASSETS.Rd +ffdacc06ed47634d2573dc65a14136dd *man/plot-binning.Rd +2db0e64a1e667bd5f26244ad866df0fb *man/plot-boxplot.Rd +e6f29eba4795f58fd09d4b01ddc39508 *man/plot-ellipses.Rd +aec6c13e84a1096d08220d425c5eec1d *man/plot-hist.Rd +0903f8cd9b38c5c7ba56824d4131cd74 *man/plot-mst.Rd +805d632b80e2e8f8da8c5c3d94dfc431 *man/plot-pairs.Rd +c4f265a109ff4657f6ec59726c44366e *man/plot-qqplot.Rd +db15ac0a3bb9c6f433286b7593d4e679 *man/plot-risk.Rd +e4d38743d705c07b14e84d10d0bbc64b *man/plot-series.Rd +1148425a4430f8c3ea914a18d88f9f5c *man/plot-similarity.Rd +42e3c47f85165de162d9bd624f36e6ad *man/plot-stars.Rd +f407709493ecda7ebf27420ef6bc556a *src/Makevars +6f7a409f4219ce1394bf16c7ce3a0f2f *src/ecodist.c +b4f61bb5bcfa23efe41ce5f84a3218fb *src/energy.c +6b26e5955f3e973aed8c110a80aee475 *tests/doRUnit.R diff -Nru fassets-2100.78/NAMESPACE fassets-2110.79/NAMESPACE --- fassets-2100.78/NAMESPACE 2009-09-24 19:08:04.000000000 +0000 +++ fassets-2110.79/NAMESPACE 2012-09-24 07:27:08.000000000 +0000 @@ -37,89 +37,67 @@ ################################################ export( + ".DEoptim", + ".MCDMeanCov", + ".OGKMeanCov", ".abcArrange", ".arwMeanCov", - "assetsArrange", - "assetsBasicStatsPlot", - "assetsBoxPercentilePlot", - "assetsBoxPlot", - "assetsBoxStatsPlot", - "assetsCorEigenPlot", - "assetsCorgramPlot", - "assetsCorImagePlot", - "assetsCorTestPlot", - "assetsCumulatedPlot", - "assetsDendrogramPlot", - "assetsFit", - "assetsHistPairsPlot", - "assetsHistPlot", - "assetsLogDensityPlot", - "assetsLPM", - "assetsMeanCov", - "assetsMomentsPlot", - "assetsNIGFitPlot", - "assetsNIGShapeTrianglePlot", - "assetsOutliers", - "assetsPairsPlot", - "assetsQQNormPlot", - "assetsReturnPlot", - "assetsRiskReturnPlot", - "assetsSelect", - "assetsSeriesPlot", - "assetsSim", - "assetsStarsPlot", ".assetsStats", - "assetsTest", - "assetsTreePlot", ".bag.fun", ".baggedMeanCov", ".bayesSteinMeanCov", + ".binaryDist", + ".braycurtisDist", + ".canberraDist", ".col.corrgram", ".cor.bagged", ".cor.mean.tawny", - ".corrgram", ".cor.shrink", + ".corDist", + ".corrgram", ".cortestPanel", ".cov.arw", ".cov.bagged", ".cov.donostah", - "covEllipsesPlot", - ".covMeanCov", - ".cov.nne.dDk", ".cov.nne.Mtovec", + ".cov.nne.dDk", ".cov.nne.nclean.sub", ".cov.nne.splusNN", ".cov.nne.vectoM", ".cov.nnve", ".cov.prior.cc", ".cov.prior.identity", - ".covRob.control", ".cov.sample.tawny", ".cov.shrink", ".cov.shrink.tawny", + ".covMeanCov", + ".covRob.control", ".denoise", - ".DEoptim", ".deoptimPlot", ".deoptimSummary", + ".differenceDist", ".dmp", ".donostahMeanCov", ".dutchPortfolioData", + ".ecodist", ".ellipsePanel", + ".euclideanDist", ".filter.RMT", - "getCenterRob", ".getCorFilter.Shrinkage", - "getCovRob", ".hclustArrange", ".hclustSelect", - ".hist", ".histPanel", + ".jaccardDist", + ".kendallDist", ".kmeansSelect", - "lambdaCVaR", ".ledoitWolfMeanCov", ".loglevel.tawny", ".lowessPanel", + ".mahalanobisDist", + ".manhattanDist", + ".maximumDist", ".mcdMeanCov", - ".MCDMeanCov", + ".minkowskiDist", ".minmaxPanel", ".mp.density.kernel", ".mp.fit.kernel", @@ -128,17 +106,16 @@ ".mp.theory", ".mst", ".mstPlot", + ".mutinfoDist", ".mveMeanCov", ".mvenergyTest", ".mvnorm.e", ".mvnormFit", ".mvshapiroTest", - "mvsnormFit", ".mvstFit", ".nnveMeanCov", ".nsca", ".numberPanel", - ".OGKMeanCov", ".orderArrange", ".panel.copula", ".panel.ellipse", @@ -149,14 +126,6 @@ ".panel.shade", ".panel.txt", ".pcaArrange", - "pfolioCVaR", - "pfolioCVaRplus", - "pfolioHist", - "pfolioMaxLoss", - "pfolioReturn", - "pfolioTargetReturn", - "pfolioTargetRisk", - "pfolioVaR", ".piePanel", ".piePtsPanel", ".ptsPanel", @@ -164,19 +133,70 @@ ".robust.cov.boot", ".sampleArrange", ".shadePanel", + ".shrinkMeanCov", ".shrinkage.c", ".shrinkage.intensity", ".shrinkage.p", ".shrinkage.r", - ".shrinkMeanCov", ".sm132PortfolioData", ".sm2vec", ".smindexes", + ".sorensenDist", ".sortIndexMST", + ".spearmanDist", ".statsArrange", ".studentMeanCov", ".txtPanel", ".usPortfolioData", ".varcov", ".vec2sm", - ".worldIndexData" ) + ".worldIndexData", + "assetsArrange", + "assetsBasicStatsPlot", + "assetsBoxPercentilePlot", + "assetsBoxPlot", + "assetsBoxStatsPlot", + "assetsCorEigenPlot", + "assetsCorImagePlot", + "assetsCorTestPlot", + "assetsCorgramPlot", + "assetsCumulatedPlot", + "assetsDendrogramPlot", + "assetsFit", + "assetsHistPairsPlot", + "assetsHistPlot", + "assetsLPM", + "assetsLogDensityPlot", + "assetsMeanCov", + "assetsMomentsPlot", + "assetsNIGFitPlot", + "assetsNIGShapeTrianglePlot", + "assetsOutliers", + "assetsPairsPlot", + "assetsQQNormPlot", + "assetsReturnPlot", + "assetsRiskReturnPlot", + "assetsSelect", + "assetsSeriesPlot", + "assetsSim", + "assetsStarsPlot", + "assetsTest", + "assetsTreePlot", + "covEllipsesPlot", + "covarRisk", + "getCenterRob", + "getCovRob", + "lambdaCVaR", + "mcr", + "mcrBeta", + "mvsnormFit", + "pfolioCVaR", + "pfolioCVaRplus", + "pfolioHist", + "pfolioMaxLoss", + "pfolioReturn", + "pfolioTargetReturn", + "pfolioTargetRisk", + "pfolioVaR", + "riskBudgets", + "riskContributions" ) diff -Nru fassets-2100.78/R/assets-arrange.R fassets-2110.79/R/assets-arrange.R --- fassets-2100.78/R/assets-arrange.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/assets-arrange.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,233 @@ + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# assetsArrange Rearranges the columns in a deta set of assets +# .pcaArrange Returns PCA correlation ordered column names +# .hclustArrange Returns hierarchical clustered column names +# .abcArrange Returns sorted column names +# .orderArrange Returns ordered column names +# .sampleArrange Returns sampled column names +# .statsArrange Returns statistically rearranged column names +################################################################################ + + +assetsArrange <- + function(x, method = c("pca", "hclust", "abc"), ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Returns ordered column names of a time Series + + # Arguments: + # x - S4 object of class 'timeSeries' + + # FUNCTION: + + # Settings: + method = match.arg(method) + FUN = paste(".", method, "Arrange", sep = "") + arrange = match.fun(FUN) + + # Return Value: + arrange(x, ...) +} + +# ------------------------------------------------------------------------------ + + +.pcaArrange <- + function(x, robust = FALSE, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Returns PCA correlation ordered column names + + # Arguments: + # x - S4 object of class 'timeSeries' + + # FUNCTION: + + # Order: + if (robust) { + if (require(robustbase)) + x.cor = covMcd(as.matrix(x), cor = TRUE, ...)$cor + else + stop("package \"robustbase\" cannot be loaded") + } else { + x.cor = cor(as.matrix(x), ...) + } + x.eigen = eigen(x.cor)$vectors[,1:2] + e1 = x.eigen[, 1] + e2 = x.eigen[, 2] + Order = order(ifelse(e1 > 0, atan(e2/e1), atan(e2/e1)+pi)) + ans = colnames(as.matrix(x))[Order] + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.hclustArrange <- + function(x, method = c("euclidean", "complete"), ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Returns hierarchical clustered column names + + # Arguments: + # x - S4 object of class 'timeSeries' + # ... + # method - the agglomeration method to be used. This should + # be (an unambiguous abbreviation of) one of "ward", "single", + # "complete", "average", "mcquitty", "median" or "centroid". + + # FUNCTION: + + # Order: + Order = hclust(dist(t(as.matrix(x)), + method = method[1]), method = method[2], ...)$order + ans = colnames(as.matrix(x))[Order] + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.abcArrange <- + function(x, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Returns sorted column names of a time Series + + # Arguments: + # x - S4 object of class 'timeSeries' + + # FUNCTION: + + # Sort: + ans = sort(colnames(as.matrix(x)), ...) + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.orderArrange <- + function(x, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Returns ordered column names of a time Series + + # Arguments: + # x - S4 object of class 'timeSeries' + + # FUNCTION: + + # Order: + ans = order(colnames(as.matrix(x)), ...) + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.sampleArrange <- + function(x, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Returns sampled column names of a time Series + + # Arguments: + # x - S4 object of class 'timeSeries' + + # FUNCTION: + + # Sample: + ans = sample(colnames(as.matrix(x)), ...) + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.statsArrange <- + function(x, FUN = colMeans, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Returns statistically rearranged column names + + # Arguments: + # x - S4 object of class 'timeSeries' + + # Note: + # Example of function Candidates: + # colStats calculates column statistics, + # colSums calculates column sums, + # colMeans calculates column means, + # colSds calculates column standard deviations, + # colVars calculates column variances, + # colSkewness calculates column skewness, + # colKurtosis calculates column kurtosis, + # colMaxs calculates maximum values in each column, + # colMins calculates minimum values in each column, + # colProds computes product of all values in each column, + # colQuantiles computes quantiles of each column. + + # FUNCTION: + + # Apply colStats Function: + fun = match.fun(FUN) + Sort = sort(fun(x, ...)) + Order = names(Sort) + ans = colnames(as.matrix(x)[, Order]) + attr(ans, "control") <- Sort + + # Return Value: + ans +} + + +################################################################################ + diff -Nru fassets-2100.78/R/assets-fit.R fassets-2110.79/R/assets-fit.R --- fassets-2100.78/R/assets-fit.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/assets-fit.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,363 @@ + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# assetsFit Fits the parameters of a set of assets +# .mvnormFit Fits a multivariate Normal distribution +# .mvsnormFit Fits a multivariate skew-Normal distribution +# .mvstFit Fits a multivariate skew-Student-t distribution +################################################################################ + + +assetsFit = + function(x, method = c("st", "snorm", "norm"), title = NULL, + description = NULL, fixed.df = NA, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Fits the parameters of a multivariate data set of assets + # and returns a list with the values for the mean, the covariance, + # the skewness, and the fatness of tails. + + # Arguments: + # x - A multivariate time series, a data frame, or any other + # rectangular object of assets which can be converted into + # a matrix by the function as.matrix. Optional Dates are + # rownames, instrument names are column names. + # type - Which type of distribution should be fitted? + # a) norm - multivariate Normal + # b) snorm - multivariate skew-Normal + # c) st - multivariate skew-Student-t + + # Value: + # The function returns a list with the following entries: + # mu - Mean values of each asset time series + # Omega - Covariance matrix of assets + # alpha - Skewness vector + # df - Degrees of freedom, measures kurtosis + + # Notes: + # Requires function "msn.mle" ans "mst.mle" from R's GPL licensed + # contributed package "sn", (C) 1998-2004 A. Azzalini. + # The list returned by this function can serve as input for the + # function assetsSim(). + + # FUNCTION: + + # Settings: + assets = as.matrix(x) + method = method[1] + colNames = colnames(x) + + # Normal Distribution: + if (method == "norm") { + # Fit Normal: + fit = list() + mu = apply(assets, 2, mean) + Omega = cov(assets) + alpha = rep(0, times = length(mu)) + df = Inf + } + + # Skew-Normal Distribution: + if (method == "snorm") { + # Fit skew-Normal: + fit = mvFit(assets, method = "snorm", ...) + mu = as.vector(fit@fit$dp$beta) + Omega = fit@fit$dp$Omega + alpha = as.vector(fit@fit$dp$alpha) + df = Inf + fit = fit@fit + } + + # Skew-Student-t Distribution: + if (method == "st") { + # Fit skew-Student: + fit = mvFit(assets, method = "st", fixed.df = fixed.df, ...) + mu = as.vector(fit@fit$beta) + Omega = fit@fit$dp$Omega + alpha = as.vector(fit@fit$dp$alpha) + df = fit@fit$dp$df + fit = fit@fit + } + + # Add Names: + names(mu) = colNames + names(alpha) = colNames + rownames(Omega) = colNames + colnames(Omega) = colNames + + # Add Title: + if (is.null(title)) + title = paste("Fitted Asset Data Model: ", method) + + # Add Description: + if (is.null(description)) + description = description() + + # Return Value: + new("fASSETS", + call = as.call(match.call()), + method = as.character(method), + model = list(mu = mu, Omega = Omega, alpha = alpha, df = df), + data = as.data.frame(x), + fit = as.list(fit), + title = as.character(title), + description = as.character(description) + ) +} + + +################################################################################ + + +.mvnormFit = + function(x, title = NULL, description = NULL, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Fits a multivariate Normal distribution + + # Arguments: + # x - A multivariate time series, a data frame, or any other + # rectangular object of assets which can be converted into + # a matrix by the function as.matrix. Optional Dates are + # rownames, instrument names are column names. + # type - Which type of distribution should be fitted? + # a) norm - multivariate Normal + # b) snorm - multivariate skew-Normal + # c) st - multivariate skew-Student-t + + # Value: + # The function returns a list with the following entries: + # mu - Mean values of each asset time series + # Omega - Covariance matrix of assets + + # Notes: + # Requires function "msn.mle" ans "mst.mle" from R's GPL licensed + # contributed package "sn", (C) 1998-2004 A. Azzalini. + # The list returned by this function can serve as input for the + # function assetsSim(). + + # FUNCTION: + + # Settings: + assets = as.matrix(x) + method = method[1] + colNames = colnames(x) + + # Fit mvNormal: + fit = list() + mu = apply(assets, 2, mean) + Omega = cov(assets) + alpha = rep(0, times = length(mu)) + df = Inf + + # Add Names: + names(mu) = colNames + names(alpha) = colNames + rownames(Omega) = colNames + colnames(Omega) = colNames + + # Add Title: + if (is.null(title)) + title = paste("Fitted Asset Data Model: ", method) + + # Add Description: + if (is.null(description)) + description = description() + + # Return Value: + new("fASSETS", + call = as.call(match.call()), + method = as.character(method), + model = list(mu = mu, Omega = Omega, alpha = alpha, df = df), + data = as.data.frame(x), + fit = as.list(fit), + title = as.character(title), + description = as.character(description) + ) +} + + +# ------------------------------------------------------------------------------ + + +mvsnormFit = + function(x, title = NULL, description = NULL, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Fits a multivariate skew-Normal distribution + + # Arguments: + # x - A multivariate time series, a data frame, or any other + # rectangular object of assets which can be converted into + # a matrix by the function as.matrix. Optional Dates are + # rownames, instrument names are column names. + # type - Which type of distribution should be fitted? + # a) norm - multivariate Normal + # b) snorm - multivariate skew-Normal + # c) st - multivariate skew-Student-t + + # Value: + # The function returns a list with the following entries: + # mu - Mean values of each asset time series + # Omega - Covariance matrix of assets + # alpha - Skewness vector + + # Notes: + # Requires function "msn.mle" ans "mst.mle" from R's GPL licensed + # contributed package "sn", (C) 1998-2004 A. Azzalini. + # The list returned by this function can serve as input for the + # function assetsSim(). + + # FUNCTION: + + # Settings: + assets = as.matrix(x) + method = method[1] + colNames = colnames(x) + + # Fit skew-Normal: + fit = mvFit(assets, method = "snorm", ...) + mu = as.vector(fit@fit$dp$beta) + Omega = fit@fit$dp$Omega + alpha = as.vector(fit@fit$dp$alpha) + df = Inf + fit = fit@fit + + + # Add Names: + names(mu) = colNames + names(alpha) = colNames + rownames(Omega) = colNames + colnames(Omega) = colNames + + # Add Title: + if (is.null(title)) + title = paste("Fitted Asset Data Model: ", method) + + # Add Description: + if (is.null(description)) + description = description() + + # Return Value: + new("fASSETS", + call = as.call(match.call()), + method = as.character(method), + model = list(mu = mu, Omega = Omega, alpha = alpha, df = df), + data = as.data.frame(x), + fit = as.list(fit), + title = as.character(title), + description = as.character(description) + ) +} + + +# ------------------------------------------------------------------------------ + + +.mvstFit = + function(x, title = NULL, description = NULL, fixed.df = NA, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Fits a multivariate skew-Student-t distribution + + # Arguments: + # x - A multivariate time series, a data frame, or any other + # rectangular object of assets which can be converted into + # a matrix by the function as.matrix. Optional Dates are + # rownames, instrument names are column names. + # type - Which type of distribution should be fitted? + # a) norm - multivariate Normal + # b) snorm - multivariate skew-Normal + # c) st - multivariate skew-Student-t + + # Value: + # The function returns a list with the following entries: + # mu - Mean values of each asset time series + # Omega - Covariance matrix of assets + # alpha - Skewness vector + # df - Degrees of freedom, measures kurtosis + + # Notes: + # Requires function "msn.mle" ans "mst.mle" from R's GPL licensed + # contributed package "sn", (C) 1998-2004 A. Azzalini. + # The list returned by this function can serve as input for the + # function assetsSim(). + + # FUNCTION: + + # Settings: + assets = as.matrix(x) + method = method[1] + colNames = colnames(x) + + # Fit skew-Student: + fit = mvFit(assets, method = "st", fixed.df = fixed.df, ...) + mu = as.vector(fit@fit$beta) + Omega = fit@fit$dp$Omega + alpha = as.vector(fit@fit$dp$alpha) + df = fit@fit$dp$df + fit = fit@fit + + # Add Names: + names(mu) = colNames + names(alpha) = colNames + rownames(Omega) = colNames + colnames(Omega) = colNames + + # Add Title: + if (is.null(title)) + title = paste("Fitted Asset Data Model: ", method) + + # Add Description: + if (is.null(description)) + description = description() + + # Return Value: + new("fASSETS", + call = as.call(match.call()), + method = as.character(method), + model = list(mu = mu, Omega = Omega, alpha = alpha, df = df), + data = as.data.frame(x), + fit = as.list(fit), + title = as.character(title), + description = as.character(description) + ) +} + + +################################################################################ + + + + + + + + + + + diff -Nru fassets-2100.78/R/assets-lpm.R fassets-2110.79/R/assets-lpm.R --- fassets-2100.78/R/assets-lpm.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/assets-lpm.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,90 @@ + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# assetsLPM Computes Asymmetric Lower Partial Moments +################################################################################ + + +assetsLPM = + function(x, tau = colMeans(x), a = 1.5, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes LPM and CLPM from multivariate time series + + # Arguments: + # x - a multivariate time series, a data frame, or any other + # rectangular object of assets which can be converted into + # a matrix by the function 'as.matrix'. Optional Dates are + # rownames, instrument names are column names. + + # Note: + # The output of this function can be used for portfolio + # optimization. + + # Example: + # LPP = as.timeSeries(data(LPP2005REC))[, 1:6] + # assetsLPM(LPP) + + # FUNCTION: + + # Transform Input: + x.mat = as.matrix(x) + nCol = ncol(x) + nRow = nrow(x) + Tau = matrix(rep(tau, nRow), byrow = TRUE, ncol = nCol) + TauX = Tau-x + X.max = ((TauX) + abs(TauX))/2 + + # Compute Lower Partial Moments: + LPM = colMeans(X.max^a) + + # Compute co-LPMs: + CLPM = diag(0, nCol) + if (a > 1) { + for (i in 1:nCol) { + for (j in 1:nCol) { + CLPM[i, j] = + mean( (X.max[, i])^(a-1) * TauX[, j] ) + } + CLPM[i, i] = LPM[i] + } + } else if (a == 1) { + for (i in 1:nCol) { + for (j in 1:nCol) { + CLPM[i, j] = + mean( sign( X.max[, i]) * TauX[, j] ) + } + CLPM[i, i] = LPM[i] + } + } + + # Result: + ans = list(mu = LPM, Sigma = CLPM) + attr(ans, "control") <- c(a = a, tau = tau) + + # Return Value: + ans +} + + +################################################################################ + + \ No newline at end of file diff -Nru fassets-2100.78/R/assets-mcr.R fassets-2110.79/R/assets-mcr.R --- fassets-2100.78/R/assets-mcr.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/assets-mcr.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,211 @@ + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# covarRisk Computes covariance portfolio risk +# mcr Computes marginal contribution to covariance risk +# mcrBeta Computes beta, the rescaled mcr to covariance risk +# riskContributions Computes covariance risk contributions +# riskBudgets Computes covariance risk budgets +################################################################################ + + +covarRisk <- +function(data, weights = NULL, FUN = "cov", ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes covariance portfolio risk + + # Arguments: + # data - a multivariate timeSeries object of financial returns + # weights - numeric vector of portfolio weights + # FUN - a covariance estimator, which returns a matrix of + # covariance estimates, by default the sample covariance + # ... - Optional arguments passed to the function FUN + + # Example: + # .covarRisk(data) + + # FUNCTION: + + # Covariance Risk: + covFun = match.fun(FUN) + COV = covFun(data) + + # Portfolio Weights: + N = ncol(COV) + if (is.null(weights)) weights = rep(1/N, N) + names(weights) = colnames(COV) + + # Covariance Portfolio Risk: + covarRisk <- sqrt( t(weights) %*% COV %*% weights )[[1, 1]] + + # Return Value: + covarRisk +} + + +# ------------------------------------------------------------------------------ + + +mcr <- +function(data, weights = NULL, FUN = "cov", ...) +{ + # A function implemented by Diethelm Wuertz + + # Description + # Computes marginal contribution to covariance risk + + # Arguments: + # data - a multivariate timeSeries object of financial returns + # weights - numeric vector of portfolio weights + # FUN - a covariance estimator, which returns a matrix of + # covariance estimates, by default the sample covariance + # ... - Optional arguments passed to the function FUN + + # Details + # The formula are implemented according to Goldberg et al., + # see also R script assetsPfolio.R + + # References: + # Lisa Goldberg et al., Extreme Risk Management, 2009 + # Scherer and Martin, Introduction to modern portfolio Optimimization + + # Example: + # data = assetsSim(100, 6); mcr(data) + + # FUNCTION: + + # Covariance Risk: + covFun = match.fun(FUN) + COV = covFun(data) + N = ncol(data) + if (is.null(weights)) weights = rep(1/N, N) + + # Marginal Contribution to Risk + mcr <- (COV %*% weights)[, 1] / covarRisk(data, weights, FUN, ...) + names(mcr) <- colnames(data) + + # Return Value: + mcr +} + + +# ------------------------------------------------------------------------------ + + +mcrBeta <- +function(data, weights = NULL, FUN = "cov", ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes beta, the rescaled marginal contribution to covariance risk + + # Arguments: + # data - a multivariate timeSeries object of financial returns + # weights - numeric vector of portfolio weights + # FUN - a covariance estimator, which returns a matrix of + # covariance estimates, by default the sample covariance + # ... - Optional arguments passed to the function FUN + + # Example: + # mcrBeta(data) + + # FUNCTION: + + # Portfolio Beta: + beta <- mcr(data, weights, FUN = FUN, ...) / + covarRisk(data, weights, FUN = FUN, ...) + + # Return Value: + beta +} + + +# ------------------------------------------------------------------------------ + + +riskContributions <- +function(data, weights = NULL, FUN = "cov", ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes covariance risk contributions + + # Arguments: + # data - a multivariate timeSeries object of financial returns + # weights - numeric vector of portfolio weights + # FUN - a covariance estimator, which returns a matrix of + # covariance estimates, by default the sample covariance + # ... - Optional arguments passed to the function FUN + + # Example: + # riskContributions(data) + + # FUNCTION: + + # Risk Contributions: + if (is.null(weights)) { + N = ncol(data) + weights = rep(1/N, times = N) + } + riskContributions <- weights * mcr(data, weights, FUN, ...) + + # Return Value: + riskContributions +} + + +# ------------------------------------------------------------------------------ + + +riskBudgets <- +function(data, weights = NULL, FUN = "cov", ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes covariance risk budgets + + # Arguments: + # data - a multivariate timeSeries object of financial returns + # weights - numeric vector of portfolio weights + # FUN - a covariance estimator, which returns a matrix of + # covariance estimates, by default the sample covariance + # ... - Optional arguments passed to the function FUN + + # Example: + # data = 100*LPP2005.RET[, 1:6]; riskBudgets(data) + + # FUNCTION: + + # Risk Budgets: + riskBudgets <- riskContributions(data, weights, FUN, ...) / + covarRisk(data, weights, FUN, ...) + + # Return Value: + riskBudgets +} + + +################################################################################ + diff -Nru fassets-2100.78/R/assets-meancov.R fassets-2110.79/R/assets-meancov.R --- fassets-2100.78/R/assets-meancov.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/assets-meancov.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,682 @@ + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# assetsMeanCov Estimates mean and variance for a set of assets +# FUNCTION: DESCRIPTION: +# .covMeanCov uses sample covariance estimation +# .mveMeanCov uses "cov.mve" from [MASS] +# .mcdMeanCov uses "cov.mcd" from [MASS] +# .studentMeanCov uses "cov.trob" from [MASS] +# .MCDMeanCov requires "covMcd" from [robustbase] +# .OGKMeanCov requires "covOGK" from [robustbase] +# .nnveMeanCov uses builtin from [covRobust] +# .shrinkMeanCov uses builtin from [corpcor] +# .baggedMeanCov uses builtin from [corpcor] +# .arwMeanCov uses builtin from [mvoutlier] +# .donostahMeanCov uses builtin from [robust] +# .bayesSteinMeanCov copy from Alexios Ghalanos +# .ledoitWolfMeanCov uses builtin from [tawny] +# .rmtMeanCov uses builtin from [tawny] +# FUNCTION: DESCRIPTION: +# getCenterRob Extracts the robust estimate for the center +# getCovRob Extracts the robust estimate for the covariance +################################################################################ + + +assetsMeanCov <- +function(x, + method = c("cov", "mve", "mcd", "MCD", "OGK", "nnve", "shrink", "bagged"), + check = TRUE, force = TRUE, baggedR = 100, sigmamu = scaleTau2, alpha = 1/2, + ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes robust mean and covariance from multivariate time series + + # Arguments: + # x - a multivariate time series, a data frame, or any other + # rectangular object of assets which can be converted into + # a matrix by the function 'as.matrix'. Optional Dates are + # rownames, instrument names are column names. + # method - Which method should be used to compute the covarinace? + # method = "cov" sample covariance computation + # method = "mve" uses "mve" from [MASS] + # method = "mcd" uses "mcd" from [MASS] + # method = "MCD" uses "MCD" from [robustbase] + # method = "OGK" uses "OGK" from [robustbase] + # method = "nnve" uses "nnve" from [covRobust] + # method = "shrink" uses "shrinkage" from [corpcor] + # method = "bagged" uses "bagging" [corpcor] + # alpha - MCD: numeric parameter controlling the size of the subsets + # over which the determinant is minimized, i.e., alpha*n observations + # are used for computing the determinant. Allowed values are between + # 0.5 and 1 and the default is 0.5. + # sigma.mu - OGK: a function that computes univariate robust location + # and scale estimates. By default it should return a single numeric + # value containing the robust scale (standard deviation) estimate. + # When mu.too is true, sigmamu() should return a numeric vector of + # length 2 containing robust location and scale estimates. See + # scaleTau2, s_Qn, s_Sn, s_mad or s_IQR for examples to be used as + # sigmamu argument. + + # Note: + # The output of this function can be used for portfolio + # optimization. + + # Example: + # DJ = 100 * returns(as.timeSeries(data(DowJones30))) + # DJ = DJ[, c("CAT", "IBM", "GE", "JPM")] + # Sample Covariance: + # assetsMeanCov(DJ, "cov") + # MASS: + # assetsMeanCov(DJ, "mve") + # assetsMeanCov(DJ, "mcd") + # require(robustbase) + # assetsMeanCov(DJ, "MCD") + # assetsMeanCov(DJ, "OGK") + # require(covRobust) + # assetsMeanCov(DJ, "nnve") + + # FUNCTION: + + # Transform Input: + x.mat = as.matrix(x) + + # Do not use: method = match.arg(method) + method = method[1] + N = ncol(x) + assetNames = colnames(x.mat) + + # Attribute Control List: + control = c(method = method[1]) + user = TRUE + + # Compute Classical Covariance: + if (method == "cov") { + # Classical Covariance Estimation: + ans = list(center = colMeans(x.mat), cov = cov(x.mat)) + user = FALSE + } + + # From R Package "robustbase": + if (method == "MCD" | method == "Mcd") { + ans = robustbase::covMcd(x.mat, alpha = alpha, ...) + mu = ans$center + Sigma = ans$cov + user = FALSE + } + if (method == "OGK" | method == "Ogk") { + ans = robustbase::covOGK(x.mat, sigmamu = scaleTau2, ...) + user = FALSE + } + + # [MASS] mve and mcd Routines: + if (method == "mve") { + # require(MASS) + ans = MASS::cov.rob(x = x.mat, method = "mve") + user = FALSE + } + if (method == "mcd") { + # require(MASS) + ans = MASS::cov.rob(x = x.mat, method = "mcd") + user = FALSE + } + + # [corpcor] Shrinkage and Bagging Routines + if (method == "shrink") { + fit = .cov.shrink(x = x.mat, ...) + ans = list(center = colMeans(x.mat), cov = fit) + user = FALSE + } + if (method == "bagged") { + fit = .cov.bagged(x = x.mat, R = baggedR, ...) + ans = list(center = colMeans(x.mat), cov = fit) + control = c(control, R = as.character(baggedR)) + user = FALSE + } + + # Nearest Neighbour Variance Estimation: + if (method == "nnve") { + fit = .cov.nnve(datamat = x.mat, ...) + ans = list(center = colMeans(x.mat), cov = fit$cov) + user = FALSE + } + + # User specified estimator: + if(user) { + fun = match.fun(method[1]) + ans = fun(x.mat, ...) + } + + # Result: + mu = center = ans$center + Sigma = cov = ans$cov + + # Add Size to Control List: + control = c(control, size = as.character(N)) + + # Add Names for Covariance Matrix to Control List: + names(mu) = assetNames + colnames(Sigma) = rownames(Sigma) = colNames = assetNames + + # Check Positive Definiteness: + if (check) { + result = isPositiveDefinite(Sigma) + if(result) { + control = c(control, posdef = "TRUE") + } else { + control = c(control, posdef = "FALSE") + } + } + + # Check Positive Definiteness: + control = c(control, forced = "FALSE") + if (force) { + control = c(control, forced = "TRUE") + if (!result) Sigma = makePositiveDefinite(Sigma) + } + + # Result: + ans = list(center = mu, cov = Sigma, mu = mu, Sigma = Sigma) + attr(ans, "control") = control + + # Return Value: + ans +} + + +################################################################################ + + +.covMeanCov <- +function(x, ...) +{ + # Description: + # Uses sample covariance estimation + + # Arguments: + # x - an object of class timeSeries + + # FUNCTION: + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + + ans = list(center = colMeans(x.mat), cov = cov(x.mat)) + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.mveMeanCov <- +function(x, ...) +{ + # Description: + + # Arguments: + # x - an object of class timeSeries + + # FUNCTION: + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + + ans = MASS::cov.rob(x = x.mat, method = "mve") + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.mcdMeanCov <- +function(x, ...) +{ + # Description: + + # Arguments: + # x - an object of class timeSeries + + # FUNCTION: + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + + ans = MASS::cov.rob(x = x.mat, method = "mcd") + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.studentMeanCov <- +function(x, ...) +{ + # Description: + + # Arguments: + # x - an object of class timeSeries + + # FUNCTION: + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + + ans = MASS::cov.trob(x, ...) + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.MCDMeanCov <- +function(x, alpha = 1/2, ...) +{ + # Description: + + # Arguments: + # x - an object of class timeSeries + + # FUNCTION: + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + + ans = robustbase::covMcd(x.mat, alpha = alpha, ...) + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.OGKMeanCov <- +function(x, sigmamu = scaleTau2, ...) +{ + # Description: + + # Arguments: + # x - an object of class timeSeries + + # FUNCTION: + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + + ans = robustbase::covOGK(x.mat, sigmamu = scaleTau2, ...) + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.nnveMeanCov <- +function(x, ...) +{ + # Description: + + # Arguments: + # x - an object of class timeSeries + + # FUNCTION: + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + fit = .cov.nnve(datamat = x.mat, ...) + + ans = list(center = colMeans(x.mat), cov = fit$cov) + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.shrinkMeanCov <- +function(x, ...) +{ + # Description: + + # Arguments: + # x - an object of class timeSeries + + # Note: + # Based on a function borrowed from package corpcor + + # FUNCTION: + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + fit = .cov.shrink(x = x.mat, ...) + + ans = list(center = colMeans(x.mat), cov = fit) + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.baggedMeanCov <- +function(x, baggedR = 100, ...) +{ + # Description: + + # Arguments: + # x - an object of class timeSeries + + # Note: + # Based on a function borrowed from package corpcor + + # FUNCTION: + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + fit = .cov.bagged(x = x.mat, R = baggedR, ...) + + ans = list(center = colMeans(x.mat), cov = fit) + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.arwMeanCov <- +function(x, ...) +{ + # Description: + # Adaptive reweighted estimator for multivariate location and scatter + # with hard-rejection weights and delta = chi2inv(1-d,p) + + # Arguments: + # x - an object of class timeSeries + + # Note: + # Based on a function borrowed from package mvoutlier + + # FUNCTION: + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + fit = .cov.arw(x = x.mat, center = colMeans(x.mat), cov = cov(x),, ...) + + ans = list(center = fit$center, cov = fit$cov) + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.donostahMeanCov <- +function(x, ...) +{ + # Description: + + # Arguments: + # x - an object of class timeSeries + + # Note: + # Based on a function borrowed from package robust + + # Settings: + x.mat = as.matrix(x) + N = ncol(x) + assetNames = colnames(x) + + ans = .cov.donostah(x = x.mat) + names(ans$center) = assetNames + rownames(ans$cov) = colnames(ans$cov) = assetNames + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.bayesSteinMeanCov <- +function(x, ...) +{ + # Description: + + # Arguments: + # x - an object of class timeSeries + + # Note: + # Based on a function written by Alexios Ghalanos + + # Bayes Stein estimator + # Alexios Ghalanos 2008 + # alexios at 4dscape.com + # This function encapsulates an example of shrinking the returns + # and covariance using Bayes-Stein shrinkage as described in + # Jorion, 1986. + + # Settings: + data <- getDataPart(x) + mu <- as.matrix(apply(data,2, FUN = function(x) mean(x))) + S <- cov(data) + k <- dim(data)[2] + n <- dim(data)[1] + one <- as.matrix(rep(1, k)) + a <- solve(S, one) + + # Constant non informative prior + mu.prior <- one * as.numeric(t(mu) %*% a/t(one) %*% a) + S.inv <- solve(S) + d <- t(mu-mu.prior) %*% S.inv %*% (mu-mu.prior) + d <- as.numeric(d) + lambda <- (k+2) / d + w <- lambda / (n+lambda) + mu.pred <- (1-w) * mu + w * mu.prior + + wc1 <- 1 / (n+lambda) + wc2 <- lambda*(n-1) / (n*(n+lambda)*(n-k-2)) + wc2 <- wc2 / as.numeric(t(one) %*% a) + V.post <- wc1 * S + wc2 * one %*% t(one) + V.pred <- S + V.post + sigma.post <- sqrt(diag(V.post)) + sigma.pred <- sqrt(diag(V.pred)) + + result <- list( + mu = mu, mu.prior = mu.prior, mu.predict = mu.pred, + V = S, V.post = V.post, V.pred = V.pred, Sigma = sqrt(diag(S)), + Sigma.post = sigma.post, Sigma.predict = sigma.pred) + + ans = list(center = result$mu.pred[,1], cov = result$V.pred) + names(ans$center) = colnames(x) + rownames(ans$cov) = colnames(ans$cov) = colnames(x) + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.ledoitWolfMeanCov <- +function(x, ...) +{ + # Description: + # Perform shrinkage on a sample covariance towards a biased covariance + + # Arguments: + # x - an object of class timeSeries + + # Details: + # This performs a covariance shrinkage estimation as specified in + # Ledoit and Wolf. Using within the larger framework only requires + # using the getCorFilter.Shrinkage function, which handles the work + # of constructing a shrinkage estimate of the covariance matrix of + # returns (and consequently its corresponding correlation matrix). + + # Note: + # Based on a function borrowed from package tawny + # Author: Brian Lee Yung Rowe + + # Settings: + data = getDataPart(x) + center = colMeans(data) + cov = .cov.shrink.tawny(data, ...) + + ans = list(center = center, cov = cov) + names(ans$center) = colnames(x) + rownames(ans$cov) = colnames(ans$cov) = colnames(x) + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +.rmtMeanCov <- +function(x, ...) +{ + # Description: + # Perform Random Matrix Theory on correlation matrix + + # Arguments: + # x - an object of class timeSeries + + # Author: + # tawnyBrian Lee Yung Rowe + + # Note: + # Based on a function borrowed from package tawny + # Author: Brian Lee Yung Rowe + + # FUNCTION: + + # Settings: + data = getDataPart(x) + center = colMeans(data) + cor = .filter.RMT(data, trace = FALSE, doplot = FALSE) + + g = colSds(data) + N = length(g) + cov = 0*cor + for (i in 1:N) + for (j in i:N) + cov[i,j] = cov[j,i] = g[i]*cor[i,j]*g[j] + + ans = list(center = center, cov = cov) + names(ans$center) = colnames(x) + rownames(ans$cov) = colnames(ans$cov) = colnames(x) + + # Return Value: + ans +} + + +################################################################################ + + +getCenterRob <- +function(object) +{ + # Return Value: + object$center +} + + +# ------------------------------------------------------------------------------ + + +getCovRob <- +function(object) +{ + # Return Value: + object$cov +} + + +################################################################################ + diff -Nru fassets-2100.78/R/assets-outliers.R fassets-2110.79/R/assets-outliers.R --- fassets-2100.78/R/assets-outliers.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/assets-outliers.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,96 @@ + +# This library is free software, you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation, either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY, without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR Description. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library, if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# assetsOutliers Detects outliers in multivariate assets sets +################################################################################ + + +assetsOutliers <- + function (x, center, cov, ...) +{ + # An adapted copy from contributed R package mvoutlier + + # Description: + # Detects outliers in a multivariate set of assets + + # Arguments: + + # Source: + # The code concerned with the outliers is from R package "mvoutliers" + # Moritz Gschwandtner + # Peter Filzmoser + + # References: + # P. Filzmoser, R.G. Garrett, and C. Reimann (2005). + # Multivariate Outlier Detection in Exploration Geochemistry. + # Computers & Geosciences. + + # FUNCTION: + + # Check timeSeries Input: + stopifnot(is.timeSeries(x)) + tS = x + x = as.matrix(x) + + # Critical Values: + n = nrow(x) + p = ncol(x) + if (p <= 10) pcrit = (0.240 - 0.0030 * p)/sqrt(n) + if (p > 10) pcrit = (0.252 - 0.0018 * p)/sqrt(n) + delta = qchisq(0.975, p) + + # Compute Mahalanobis Squared Distances: + d2 = mahalanobis(x, center, cov) + + # Detect Outliers: + d2ord = sort(d2) + dif = pchisq(d2ord, p) - (0.5:n)/n + i = (d2ord >= delta) & (dif > 0) + if (sum(i) == 0) alfan = 0 else alfan = max(dif[i]) + if (alfan < pcrit) alfan = 0 + if (alfan > 0) cn = max(d2ord[n-ceiling(n*alfan)], delta) else cn = Inf + w = d2 < cn + m = apply(x[w, ], 2, mean) + c1 = as.matrix(x - rep(1, n) %*% t(m)) + c = (t(c1) %*% diag(w) %*% c1)/sum(w) + + # Identify Outliers: + outliers = (1:dim(x)[1])[!w] + if (length(outliers) == 0) { + outliers = NA + } else { + names(outliers) = rownames(x)[outliers] + } + + # Compose Result: + ans = list( + center = m, + cov = c, + cor = cov2cor(c), + quantile = cn, + outliers = outliers, + series = tS[outliers, ]) + + # Return Value: + ans +} + + +################################################################################ + diff -Nru fassets-2100.78/R/assets-portfolio.R fassets-2110.79/R/assets-portfolio.R --- fassets-2100.78/R/assets-portfolio.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/assets-portfolio.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,405 @@ + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR Description. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# pfolioVaR Computes VaR for a portfolio of assets +# pfolioCVaR Computes CVaR for a portfoluio of assets +# pfolioCVaRplus Computes CVaR-Plus for a portfolio of assets +# lambdaCVaR Computes CVaR's atomic split value lambda +# FUNCTION: DESCRIPTION: +# pfolioMaxLoss Computes maximum loss for a portfolio +# pfolioReturn Computes return series for a portfolio +# pfolioTargetReturn Computes target return for a portfolio +# pfolioTargetRisk Computes target risk for a portfolio +# pfolioHist Plots a histogram of portfolio returns +################################################################################ + + +pfolioVaR <- + function(x, weights = NULL, alpha = 0.05) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Compute Value-at-Risk for a portfolio of assets + + # Arguments: + # x - a time series, data.frame or any other rectangular object + # of assets which can be written as a matrix object + # weights - a numeric vector of weights + # alpha - the confidence level + + # FUNCTION: + + # Transform: + x = as.matrix(x) + + # Compute Portfolio VaR: + if (is.null(weights)) + weights = rep(1/dim(x)[[2]], dim(x)[[2]]) + n = dim(x)[1] + x = apply(t(t(x) * weights), 1, sum) + n.alpha = max(floor(n * alpha)) + ans = as.vector(sort(x)[n.alpha]) + names(ans) = "VaR" + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +pfolioCVaRplus <- + function(x, weights = NULL, alpha = 0.05) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Compute Value-at-Risk Plus for a portfolio of assets + + # Arguments: + # x - a time series, data.frame or any other rectangular object + # of assets which can be written as a matrix object + # weights - a numeric vector of weights + # alpha - the confidence level + + # FUNCTION: + + # Transform: + x = as.matrix(x) + + # Compute Portfolio CVaRplus: + if (is.null(weights)) { + weights = rep(1/dim(x)[[2]], dim(x)[[2]]) + } + n = dim(x)[1] + x = apply(t(t(x) * weights), 1, sum) + n.alpha = max(1, floor(n * alpha)-1) + ans = as.vector(mean(sort(x)[1:n.alpha])) + names(ans) = "CVaRplus" + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +pfolioCVaR <- + function(x, weights = NULL, alpha = 0.05) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Compute Conditional Value-at-risk for a portfolio of assets + + # Arguments: + # x - a time series, data.frame or any other rectangular object + # of assets which can be written as a matrix object + # weights - a numeric vector of weights + # alpha - the confidence level + # lambda - split value + + # FUNCTION: + + # Transform: + data = as.matrix(x) + + # Input Data: + if (is.null(weights)) { + weights = rep(1/dim(data)[[2]], dim(data)[[2]]) + } + n = dim(data)[1] + Rp = apply(t(t(data)*weights), 1, sum) + + # Sort the Portfolio returns Y + sorted = sort(Rp) + + # Compute Portfolio VaR: + n.alpha = floor(n*alpha) + VaR = sorted[n.alpha] + + # Compute Portfolio CVaRplus: + n.alpha = max(1, floor(n*alpha)-1) + CVaRplus = mean(sorted[1:n.alpha]) + + # Compute Portfolio CVaR: + lambda = 1 - floor(n*alpha)/(n*alpha) + ans = as.vector(lambda*VaR + (1-lambda)*CVaRplus) + names(ans) = "CVaR" + attr(ans, "control") = c(CVaRplus = CVaRplus, lambda = lambda) + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +lambdaCVaR <- + function(n, alpha = 0.05) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes CVaR's atomic split value lambda + + # Arguments: + # n - the number of oberservations + # alpha - the confidence interval + + # FUNCTION: + + # Compute CVaR lambda: + lambda = 1 - floor(alpha * n) / (alpha * n) + names(lambda) = "lambda" + + # Return Value: + lambda +} + + +################################################################################ + + +pfolioMaxLoss <- + function(x, weights = NULL) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes maximum loss for a portfolio of assets + + # Arguments: + # x - a timeSeries, data.frame or any other rectangular object + # of assets which can be written as a matrix object + # weights - the vector of weights + # alpha - the confidence level + + # FUNCTION: + + # Transform: + x = as.matrix(x) + + # Compute MaxLoss [MinReturn]: + if (is.null(weights)) { + weights = rep(1/dim(x)[[2]], dim(x)[[2]]) + } + x = apply(t(t(x)*weights), 1, sum) + ans = min(x) + + # Return Value: + ans +} + + +# ------------------------------------------------------------------------------ + + +pfolioReturn <- + function(x, weights = NULL) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes return value of a portfolio + + # Arguments: + # x - a timeSeries, data.frame or any other rectangular object + # of assets which can be written as a matrix object + # weights - the vector of weights + + # FUNCTION: + + # Time Series Object ? + if ( class(x) == "timeSeries" ) { + TS = TRUE + data = x[, 1] + } else { + TS = FALSE + } + + # Compute Portfolio Returns: + x = as.matrix(x) + if (is.null(weights)) + weights = rep(1/dim(x)[[2]], dim(x)[[2]]) + ans = t(t(apply(t(t(x)*weights), 1, sum))) + colnames(ans) = "Portfolio" + + # Time Series Object: + if (TS) { + series(data) = ans + } else { + data = ans + } + + # Return Value: + data +} + + +# ------------------------------------------------------------------------------ + + +pfolioTargetReturn <- + function(x, weights = NULL) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes return value of a portfolio + + # Arguments: + # x - a timeSeries, data.frame or any other rectangular object + # of assets which can be written as a matrix object + # weights - the vector of weights + + # FUNCTION: + + # Transform: + x = as.matrix(x) + + # Compute Portfolio Returns: + ans = mean(pfolioReturn(x = x, weights = weights)) + + # Return Value: + names(ans) = "TargetReturn" + ans +} + + +# ------------------------------------------------------------------------------ + + +pfolioTargetRisk <- + function(x, weights = NULL) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Computes risk from covariance matrix of a portfolio + + # Arguments: + # x - a timeSeries, data.frame or any other rectangular object + # of assets which can be written as a matrix object + # weights - the vector of weights + + # FUNCTION: + + # Transform: + x = as.matrix(x) + + # Compute Portfolio Returns: + if (is.null(weights)) + weights = rep(1/dim(x)[[2]], dim(x)[[2]]) + ans = as.vector(sqrt(weights %*% cov(x) %*% weights)) + + # Return Value: + names(ans) = "TargetRisk" + ans +} + + + +# ------------------------------------------------------------------------------ + + +pfolioHist <- + function(x, weights = NULL, alpha = 0.05, range = NULL, details = TRUE, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Plots a histogram of the returns of a portfolio + + # Arguments: + # x - a timeSeries, data.frame or any other rectangular object + # of assets which can be written as a matrix object + # weights - the vector of weights + + # FUNCTION: + + # Transform: + x = as.matrix(x) + + # Suppress Warnings: + opt = options() + options(warn = -1) + + # Plot Portfolio Returns: + Returns = pfolioReturn(x = x, weights = weights) + if (is.null(range)) { + lim = 1.05 * pfolioMaxLoss(x = x, weights = weights)[[1]] + xlim = c(lim, -lim) + } else { + xlim = range + } + Histogram = hist(Returns, xlim = xlim, xlab = "Portfolio Return %", + probability = TRUE, col = "steelblue4", border = "white", ...) + r = seq(xlim[1], xlim[2], length = 201) + lines(r, dnorm(r, mean = mean(Returns), sd = sd(Returns)), ...) + + points(Returns, rep(0, length(Returns)), pch = 20, + col = "orange", cex = 1.25) + + # Add VaR, CVaRplus and MaxLoss: + V1 = pfolioVaR(x = x, weights = weights, alpha = alpha)[[1]] + abline(v = V1, col = "blue", ...) + V2 = pfolioCVaRplus(x = x, weights = weights, alpha = alpha)[[1]] + abline(v = V2, col = "red", ...) + V3 = pfolioMaxLoss(x = x, weights = weights)[[1]] + abline(v = V3, col = "green", ...) + V4 = as.vector(mean(Returns))[[1]] + V5 = as.vector(sd(Returns))[[1]] + + yt = max(density(Returns)$y) + + text(V1, yt, as.character(round(V1, 2)), cex = 0.75, col = "orange") + text(V2, yt, as.character(round(V2, 2)), cex = 0.75, col = "orange") + text(V3, yt, as.character(round(V3, 2)), cex = 0.75, col = "orange") + text(V4, yt, as.character(round(V4, 2)), cex = 0.75, col = "orange") + + yt = 0.95 * yt + text(V1, yt, "VaR", cex = 0.75, col = "orange") + text(V2, yt, "CVaR+", cex = 0.75, col = "orange") + text(V3, yt, "maxLoss", cex = 0.75, col = "orange") + text(V4, yt, "Mean", cex = 0.75, col = "orange") + + # Result: + options(opt) + ans = list(VaR = V1, VaRplus = V2, maxLoss = V3, mean = V4, sd = V5) + if (details) { + cat("\nVaR: ", V1) + cat("\nVaRplus: ", V2) + cat("\nmax Loss: ", V3) + cat("\nMean: ", V4) + cat("\nStDev: ", V5) + cat("\n") + } + + # Return Value: + invisible(ans) +} + + +################################################################################ + diff -Nru fassets-2100.78/R/assets-select.R fassets-2110.79/R/assets-select.R --- fassets-2100.78/R/assets-select.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/assets-select.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,121 @@ + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# assetsSelect Selects similar or dissimilar assets +# .hclustSelect Selects due to hierarchical clustering +# .kmeansSelect Selects due to k-means clustering +################################################################################ + + +assetsSelect = + function(x, method = c("hclust", "kmeans"), control = NULL, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Clusters a set of assets + + # Arguments: + # method - which algorithm should be used? + # hclust - Hierarchical clustering on a set of dissimilarities + # kmeans - k-means clustering on a data matrix + + # FUNCTION: + + # Selection: + # do not method = match.arg(method) to allow for user specified clustering + method = method[1] + + # Transform to matrix: + if (class(x) == "timeSeries") { + x = as.matrix(x) + } + + # Compose Function: + fun = paste(".", method, "Select", sep = "") + FUN = match.fun(fun) + + # Cluster: + ans = FUN(x, control, ...) + + # Return Value: + ans +} + + +################################################################################ + + +.hclustSelect <- + function(x, control = NULL, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Hierarchical Clustering + + # FUNCTION: + + # Method: + if (is.null(control)) + control = c(measure = "euclidean", method = "complete") + measure = control[1] + method = control[2] + + # hclust: + ans = hclust(dist(t(x), method = measure), method = method, ...) + class(ans) = c("list", "hclust") + + # Return Value: + ans +} + + +################################################################################ + + +.kmeansSelect <- + function(x, control = NULL, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # kmeans Clustering + + # Note: + # centers must be specified by the user! + + # FUNCTION: + + # Method: + if (is.null(control)) + control = c(centers = 5, algorithm = "Hartigan-Wong") + centers = as.integer(control[1]) + algorithm = control[2] + + # kmeans: + ans = kmeans(x = t(x), centers = centers, algorithm = algorithm, ...) + class(ans) = c("list", "kmeans") + + # Return Value: + ans +} + + +################################################################################ \ No newline at end of file diff -Nru fassets-2100.78/R/assets-simulate.R fassets-2110.79/R/assets-simulate.R --- fassets-2100.78/R/assets-simulate.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/assets-simulate.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,89 @@ + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# assetsSim Simulates a set of artificial assets +################################################################################ + + +assetsSim <- + function(n, dim = 2, model = list(mu = rep(0, dim), Omega = diag(dim), + alpha = rep(0, dim), df = Inf), assetNames = NULL) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Simulates a multivariate set of asset log-returns distributed + # according to a Normal, skew-Normal, or skew Student-t Distribution + # Function. + + # Arguments: + # n - the number of data records to be simulated + # dim - the dimension number, i.e. the number of assets to be simulated + # model - a list with the model parameters: + # mu - the numeric vector of mean values of each asset time series + # Omega - the covariance matrix of assets + # alpha - the skewness vector + # df - the degrees of freedom, a measures for the kurtosis + # assetNames - a string vector of asset names, by default NULL + # which creates asset names as "V1", "V2", ..., "Vd", where + # d denotes the dimension + + # Notes: + # Requires function "msn.quantities" from R's GPL licensed + # contributed package "sn", (C) 1998-2004 A. Azzalini. + # The model can also be the value returned by model slot from + # function assetsFit(). + + # Example: + # assetsSim(n=25) + # assetsSim(n=25, assetNames = c("RETURN-1", "RETURN-2") + # assetsSim(n=25, list(mu=c(0,0), Omega=diag(2), alpha=c(0,0), df=4)) + + # FUNCTION: + + # Dimensions: + d = length(model$alpha) + if ( length(model$mu) != d | any(dim(model$Omega) != c(d, d))) + stop("dimensions of arguments do not match") + + # Adapted from contributed R package "sn:rmsn" + Z = msn.quantities(model$mu, model$Omega, model$alpha) + y = matrix(rnorm(n * d), n, d) %*% chol(Z$Psi) + abs.y0 = matrix(rep(abs(rnorm(n)), d), ncol = d) + z = Z$delta * t(abs.y0) + sqrt(1 - Z$delta^2) * t(y) + + # Select: + if (model$df == Inf) { + ans = t(model$mu + Z$omega * z) + } else { + x = rchisq(n, model$df)/model$df + z = t(model$mu + Z$omega * z) + ans = t(model$mu + t(sqrt(x) * z)) + } + + # Dimnames: + dimnames(ans)[[2]] = assetNames + + # Return Value: + as.data.frame(ans) +} + + +################################################################################ + diff -Nru fassets-2100.78/R/assets-test.R fassets-2110.79/R/assets-test.R --- fassets-2100.78/R/assets-test.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/assets-test.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,191 @@ + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# assetsTest Tests for multivariate Normal Assets +# .mvshapiroTest Multivariate Shapiro Test +# .mvenergyTest Multivariate E-Statistic (Energy) Test +################################################################################ + + +assetsTest = +function(x, method = c("shapiro", "energy"), Replicates = 100, +title = NULL, description = NULL) +{ + # Description: + # Tests for multivariate Normal Assets + + # Example: + # .mvnormTest(x = assetsSim(100)) + # .mvnormTest(x = assetsSim(100), method = "e", Replicates = 99) + + # FUNCTION: + + # Test: + method = match.arg(method) + if (method == "shapiro") { + test = .mvshapiroTest(x) + } + if (method == "energy") { + test = .mvenergyTest(x, Replicates = Replicates) + } + + # Return Value: + test +} + + +# ------------------------------------------------------------------------------ + + +.mvshapiroTest = +function(x, title = NULL, description = NULL) +{ + # Description: + # Computes Shapiro's normality test for multivariate variables + + # Note: + # Reimplemented function, doesn't require the contributed R package + # mvnormtest + + # Author: + # Slawomir Jarek + # License: GPL + + # Source: + # Package: mvnormtest + # Version: 0.1-6 + # Date: 2005-04-02 + # Title: Normality test for multivariate variables + # Author: Slawomir Jarek + # Maintainer: Slawomir Jarek + # Description: Generalization of shapiro-wilk test for + # multivariate variables. + + # Example: + # .mvshapiroTest(x = assetsSim(100)) + + # FUNCTION: + + # Transform: + U = t(as.matrix(x)) + + # Test: + n = ncol(U) + if (n < 3 || n > 5000) stop("sample size must be between 3 and 5000") + rng = range(U) + rng = rng[2] - rng[1] + if (rng == 0) + stop("all `U[]' are identical") + Us = apply(U, 1, mean) + R = U-Us + M.1 = solve(R %*% t(R), tol = 1e-18) + Rmax = diag(t(R) %*% M.1 %*% R) + C = M.1 %*% R[, which.max(Rmax)] + Z = t(C) %*% U + test = shapiro.test(Z) + names(test$p.value) = "" + class(test) = "list" + + # Add title and description: + if (is.null(title)) title = "Multivariate Shapiro Test" + if (is.null(description)) description = description() + + # Return Value: + new("fHTEST", + call = match.call(), + data = list(x = x), + test = test, + title = title, + description = description) +} + + +# ------------------------------------------------------------------------------ + + +.mvenergyTest = +function(x, Replicates = 99, title = NULL, description = NULL) +{ + # Description: + # Computes E-statistics test for multivariate variables + + # Note: + # Reimplemented function, doesn't require the contributed + # R package energy, we only use the C Program here. + + # Source: + # Maria L. Rizzo and + # Gabor J. Szekely + # License: GPL 2.0 or later + + # Example: + # .mvenergyTest(x = assetsSim(100), 99) + + # FUNCTION: + + # Transform: + if (class(x) == "timeSeries") x = series(x) + x = as.matrix(x) + + # Test: + R = Replicates + + # RVs: + n <- nrow(x) + d <- ncol(x) + ran.gen = function(x, y) return(matrix(rnorm(n * d), nrow = n, ncol = d)) + + # Parametric Mini Boot: + strata = rep(1, n) + n <- nrow(x) + temp.str <- strata + strata <- tapply(1:n, as.numeric(strata)) + t0 <- .mvnorm.e(x) + lt0 <- length(t0) + t.star <- matrix(NA, sum(R), lt0) + pred.i <- NULL + for(r in 1:R) t.star[r, ] <- .mvnorm.e(ran.gen(x, NULL)) + + # Result: + test <- list( + statistic = c("E-Statistic" = t0), + p.value = 1 - mean(t.star < t0), + method = "Energy Test", + data.name = paste("x, obs ", n, ", dim ", d, ", reps ", R, sep = "")) + names(test$p.value) = "" + class(test) = "list" + + # Add: + if (is.null(title)) title = test$method + if (is.null(description)) description = description() + + # Return Value: + new("fHTEST", + call = match.call(), + data = list(x = x), + test = test, + title = title, + description = description) +} + + +################################################################################ + + + diff -Nru fassets-2100.78/R/assetsArrange.R fassets-2110.79/R/assetsArrange.R --- fassets-2100.78/R/assetsArrange.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/assetsArrange.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,233 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - - -################################################################################ -# FUNCTION: DESCRIPTION: -# assetsArrange Rearranges the columns in a deta set of assets -# .pcaArrange Returns PCA correlation ordered column names -# .hclustArrange Returns hierarchical clustered column names -# .abcArrange Returns sorted column names -# .orderArrange Returns ordered column names -# .sampleArrange Returns sampled column names -# .statsArrange Returns statistically rearranged column names -################################################################################ - - -assetsArrange <- - function(x, method = c("pca", "hclust", "abc"), ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Returns ordered column names of a time Series - - # Arguments: - # x - S4 object of class 'timeSeries' - - # FUNCTION: - - # Settings: - method = match.arg(method) - FUN = paste(".", method, "Arrange", sep = "") - arrange = match.fun(FUN) - - # Return Value: - arrange(x, ...) -} - -# ------------------------------------------------------------------------------ - - -.pcaArrange <- - function(x, robust = FALSE, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Returns PCA correlation ordered column names - - # Arguments: - # x - S4 object of class 'timeSeries' - - # FUNCTION: - - # Order: - if (robust) { - if (require(robustbase)) - x.cor = covMcd(as.matrix(x), cor = TRUE, ...)$cor - else - stop("package \"robustbase\" cannot be loaded") - } else { - x.cor = cor(as.matrix(x), ...) - } - x.eigen = eigen(x.cor)$vectors[,1:2] - e1 = x.eigen[, 1] - e2 = x.eigen[, 2] - Order = order(ifelse(e1 > 0, atan(e2/e1), atan(e2/e1)+pi)) - ans = colnames(as.matrix(x))[Order] - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.hclustArrange <- - function(x, method = c("euclidean", "complete"), ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Returns hierarchical clustered column names - - # Arguments: - # x - S4 object of class 'timeSeries' - # ... - # method - the agglomeration method to be used. This should - # be (an unambiguous abbreviation of) one of "ward", "single", - # "complete", "average", "mcquitty", "median" or "centroid". - - # FUNCTION: - - # Order: - Order = hclust(dist(t(as.matrix(x)), - method = method[1]), method = method[2], ...)$order - ans = colnames(as.matrix(x))[Order] - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.abcArrange <- - function(x, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Returns sorted column names of a time Series - - # Arguments: - # x - S4 object of class 'timeSeries' - - # FUNCTION: - - # Sort: - ans = sort(colnames(as.matrix(x)), ...) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.orderArrange <- - function(x, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Returns ordered column names of a time Series - - # Arguments: - # x - S4 object of class 'timeSeries' - - # FUNCTION: - - # Order: - ans = order(colnames(as.matrix(x)), ...) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.sampleArrange <- - function(x, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Returns sampled column names of a time Series - - # Arguments: - # x - S4 object of class 'timeSeries' - - # FUNCTION: - - # Sample: - ans = sample(colnames(as.matrix(x)), ...) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.statsArrange <- - function(x, FUN = colMeans, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Returns statistically rearranged column names - - # Arguments: - # x - S4 object of class 'timeSeries' - - # Note: - # Example of function Candidates: - # colStats calculates column statistics, - # colSums calculates column sums, - # colMeans calculates column means, - # colSds calculates column standard deviations, - # colVars calculates column variances, - # colSkewness calculates column skewness, - # colKurtosis calculates column kurtosis, - # colMaxs calculates maximum values in each column, - # colMins calculates minimum values in each column, - # colProds computes product of all values in each column, - # colQuantiles computes quantiles of each column. - - # FUNCTION: - - # Apply colStats Function: - fun = match.fun(FUN) - Sort = sort(fun(x, ...)) - Order = names(Sort) - ans = colnames(as.matrix(x)[, Order]) - attr(ans, "control") <- Sort - - # Return Value: - ans -} - - -################################################################################ - diff -Nru fassets-2100.78/R/assetsFit.R fassets-2110.79/R/assetsFit.R --- fassets-2100.78/R/assetsFit.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/assetsFit.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,363 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - - -################################################################################ -# FUNCTION: DESCRIPTION: -# assetsFit Fits the parameters of a set of assets -# .mvnormFit Fits a multivariate Normal distribution -# .mvsnormFit Fits a multivariate skew-Normal distribution -# .mvstFit Fits a multivariate skew-Student-t distribution -################################################################################ - - -assetsFit = - function(x, method = c("st", "snorm", "norm"), title = NULL, - description = NULL, fixed.df = NA, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Fits the parameters of a multivariate data set of assets - # and returns a list with the values for the mean, the covariance, - # the skewness, and the fatness of tails. - - # Arguments: - # x - A multivariate time series, a data frame, or any other - # rectangular object of assets which can be converted into - # a matrix by the function as.matrix. Optional Dates are - # rownames, instrument names are column names. - # type - Which type of distribution should be fitted? - # a) norm - multivariate Normal - # b) snorm - multivariate skew-Normal - # c) st - multivariate skew-Student-t - - # Value: - # The function returns a list with the following entries: - # mu - Mean values of each asset time series - # Omega - Covariance matrix of assets - # alpha - Skewness vector - # df - Degrees of freedom, measures kurtosis - - # Notes: - # Requires function "msn.mle" ans "mst.mle" from R's GPL licensed - # contributed package "sn", (C) 1998-2004 A. Azzalini. - # The list returned by this function can serve as input for the - # function assetsSim(). - - # FUNCTION: - - # Settings: - assets = as.matrix(x) - method = method[1] - colNames = colnames(x) - - # Normal Distribution: - if (method == "norm") { - # Fit Normal: - fit = list() - mu = apply(assets, 2, mean) - Omega = cov(assets) - alpha = rep(0, times = length(mu)) - df = Inf - } - - # Skew-Normal Distribution: - if (method == "snorm") { - # Fit skew-Normal: - fit = mvFit(assets, method = "snorm", ...) - mu = as.vector(fit@fit$dp$beta) - Omega = fit@fit$dp$Omega - alpha = as.vector(fit@fit$dp$alpha) - df = Inf - fit = fit@fit - } - - # Skew-Student-t Distribution: - if (method == "st") { - # Fit skew-Student: - fit = mvFit(assets, method = "st", fixed.df = fixed.df, ...) - mu = as.vector(fit@fit$beta) - Omega = fit@fit$dp$Omega - alpha = as.vector(fit@fit$dp$alpha) - df = fit@fit$dp$df - fit = fit@fit - } - - # Add Names: - names(mu) = colNames - names(alpha) = colNames - rownames(Omega) = colNames - colnames(Omega) = colNames - - # Add Title: - if (is.null(title)) - title = paste("Fitted Asset Data Model: ", method) - - # Add Description: - if (is.null(description)) - description = description() - - # Return Value: - new("fASSETS", - call = as.call(match.call()), - method = as.character(method), - model = list(mu = mu, Omega = Omega, alpha = alpha, df = df), - data = as.data.frame(x), - fit = as.list(fit), - title = as.character(title), - description = as.character(description) - ) -} - - -################################################################################ - - -.mvnormFit = - function(x, title = NULL, description = NULL, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Fits a multivariate Normal distribution - - # Arguments: - # x - A multivariate time series, a data frame, or any other - # rectangular object of assets which can be converted into - # a matrix by the function as.matrix. Optional Dates are - # rownames, instrument names are column names. - # type - Which type of distribution should be fitted? - # a) norm - multivariate Normal - # b) snorm - multivariate skew-Normal - # c) st - multivariate skew-Student-t - - # Value: - # The function returns a list with the following entries: - # mu - Mean values of each asset time series - # Omega - Covariance matrix of assets - - # Notes: - # Requires function "msn.mle" ans "mst.mle" from R's GPL licensed - # contributed package "sn", (C) 1998-2004 A. Azzalini. - # The list returned by this function can serve as input for the - # function assetsSim(). - - # FUNCTION: - - # Settings: - assets = as.matrix(x) - method = method[1] - colNames = colnames(x) - - # Fit mvNormal: - fit = list() - mu = apply(assets, 2, mean) - Omega = cov(assets) - alpha = rep(0, times = length(mu)) - df = Inf - - # Add Names: - names(mu) = colNames - names(alpha) = colNames - rownames(Omega) = colNames - colnames(Omega) = colNames - - # Add Title: - if (is.null(title)) - title = paste("Fitted Asset Data Model: ", method) - - # Add Description: - if (is.null(description)) - description = description() - - # Return Value: - new("fASSETS", - call = as.call(match.call()), - method = as.character(method), - model = list(mu = mu, Omega = Omega, alpha = alpha, df = df), - data = as.data.frame(x), - fit = as.list(fit), - title = as.character(title), - description = as.character(description) - ) -} - - -# ------------------------------------------------------------------------------ - - -mvsnormFit = - function(x, title = NULL, description = NULL, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Fits a multivariate skew-Normal distribution - - # Arguments: - # x - A multivariate time series, a data frame, or any other - # rectangular object of assets which can be converted into - # a matrix by the function as.matrix. Optional Dates are - # rownames, instrument names are column names. - # type - Which type of distribution should be fitted? - # a) norm - multivariate Normal - # b) snorm - multivariate skew-Normal - # c) st - multivariate skew-Student-t - - # Value: - # The function returns a list with the following entries: - # mu - Mean values of each asset time series - # Omega - Covariance matrix of assets - # alpha - Skewness vector - - # Notes: - # Requires function "msn.mle" ans "mst.mle" from R's GPL licensed - # contributed package "sn", (C) 1998-2004 A. Azzalini. - # The list returned by this function can serve as input for the - # function assetsSim(). - - # FUNCTION: - - # Settings: - assets = as.matrix(x) - method = method[1] - colNames = colnames(x) - - # Fit skew-Normal: - fit = mvFit(assets, method = "snorm", ...) - mu = as.vector(fit@fit$dp$beta) - Omega = fit@fit$dp$Omega - alpha = as.vector(fit@fit$dp$alpha) - df = Inf - fit = fit@fit - - - # Add Names: - names(mu) = colNames - names(alpha) = colNames - rownames(Omega) = colNames - colnames(Omega) = colNames - - # Add Title: - if (is.null(title)) - title = paste("Fitted Asset Data Model: ", method) - - # Add Description: - if (is.null(description)) - description = description() - - # Return Value: - new("fASSETS", - call = as.call(match.call()), - method = as.character(method), - model = list(mu = mu, Omega = Omega, alpha = alpha, df = df), - data = as.data.frame(x), - fit = as.list(fit), - title = as.character(title), - description = as.character(description) - ) -} - - -# ------------------------------------------------------------------------------ - - -.mvstFit = - function(x, title = NULL, description = NULL, fixed.df = NA, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Fits a multivariate skew-Student-t distribution - - # Arguments: - # x - A multivariate time series, a data frame, or any other - # rectangular object of assets which can be converted into - # a matrix by the function as.matrix. Optional Dates are - # rownames, instrument names are column names. - # type - Which type of distribution should be fitted? - # a) norm - multivariate Normal - # b) snorm - multivariate skew-Normal - # c) st - multivariate skew-Student-t - - # Value: - # The function returns a list with the following entries: - # mu - Mean values of each asset time series - # Omega - Covariance matrix of assets - # alpha - Skewness vector - # df - Degrees of freedom, measures kurtosis - - # Notes: - # Requires function "msn.mle" ans "mst.mle" from R's GPL licensed - # contributed package "sn", (C) 1998-2004 A. Azzalini. - # The list returned by this function can serve as input for the - # function assetsSim(). - - # FUNCTION: - - # Settings: - assets = as.matrix(x) - method = method[1] - colNames = colnames(x) - - # Fit skew-Student: - fit = mvFit(assets, method = "st", fixed.df = fixed.df, ...) - mu = as.vector(fit@fit$beta) - Omega = fit@fit$dp$Omega - alpha = as.vector(fit@fit$dp$alpha) - df = fit@fit$dp$df - fit = fit@fit - - # Add Names: - names(mu) = colNames - names(alpha) = colNames - rownames(Omega) = colNames - colnames(Omega) = colNames - - # Add Title: - if (is.null(title)) - title = paste("Fitted Asset Data Model: ", method) - - # Add Description: - if (is.null(description)) - description = description() - - # Return Value: - new("fASSETS", - call = as.call(match.call()), - method = as.character(method), - model = list(mu = mu, Omega = Omega, alpha = alpha, df = df), - data = as.data.frame(x), - fit = as.list(fit), - title = as.character(title), - description = as.character(description) - ) -} - - -################################################################################ - - - - - - - - - - - diff -Nru fassets-2100.78/R/assetsLPM.R fassets-2110.79/R/assetsLPM.R --- fassets-2100.78/R/assetsLPM.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/assetsLPM.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,90 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - - -################################################################################ -# FUNCTION: DESCRIPTION: -# assetsLPM Computes Asymmetric Lower Partial Moments -################################################################################ - - -assetsLPM = - function(x, tau = colMeans(x), a = 1.5, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Computes LPM and CLPM from multivariate time series - - # Arguments: - # x - a multivariate time series, a data frame, or any other - # rectangular object of assets which can be converted into - # a matrix by the function 'as.matrix'. Optional Dates are - # rownames, instrument names are column names. - - # Note: - # The output of this function can be used for portfolio - # optimization. - - # Example: - # LPP = as.timeSeries(data(LPP2005REC))[, 1:6] - # assetsLPM(LPP) - - # FUNCTION: - - # Transform Input: - x.mat = as.matrix(x) - nCol = ncol(x) - nRow = nrow(x) - Tau = matrix(rep(tau, nRow), byrow = TRUE, ncol = nCol) - TauX = Tau-x - X.max = ((TauX) + abs(TauX))/2 - - # Compute Lower Partial Moments: - LPM = colMeans(X.max^a) - - # Compute co-LPMs: - CLPM = diag(0, nCol) - if (a > 1) { - for (i in 1:nCol) { - for (j in 1:nCol) { - CLPM[i, j] = - mean( (X.max[, i])^(a-1) * TauX[, j] ) - } - CLPM[i, i] = LPM[i] - } - } else if (a == 1) { - for (i in 1:nCol) { - for (j in 1:nCol) { - CLPM[i, j] = - mean( sign( X.max[, i]) * TauX[, j] ) - } - CLPM[i, i] = LPM[i] - } - } - - # Result: - ans = list(mu = LPM, Sigma = CLPM) - attr(ans, "control") <- c(a = a, tau = tau) - - # Return Value: - ans -} - - -################################################################################ - - \ No newline at end of file diff -Nru fassets-2100.78/R/assetsMeanCov.R fassets-2110.79/R/assetsMeanCov.R --- fassets-2100.78/R/assetsMeanCov.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/assetsMeanCov.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,682 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - - -################################################################################ -# FUNCTION: DESCRIPTION: -# assetsMeanCov Estimates mean and variance for a set of assets -# FUNCTION: DESCRIPTION: -# .covMeanCov uses sample covariance estimation -# .mveMeanCov uses "cov.mve" from [MASS] -# .mcdMeanCov uses "cov.mcd" from [MASS] -# .studentMeanCov uses "cov.trob" from [MASS] -# .MCDMeanCov requires "covMcd" from [robustbase] -# .OGKMeanCov requires "covOGK" from [robustbase] -# .nnveMeanCov uses builtin from [covRobust] -# .shrinkMeanCov uses builtin from [corpcor] -# .baggedMeanCov uses builtin from [corpcor] -# .arwMeanCov uses builtin from [mvoutlier] -# .donostahMeanCov uses builtin from [robust] -# .bayesSteinMeanCov copy from Alexios Ghalanos -# .ledoitWolfMeanCov uses builtin from [tawny] -# .rmtMeanCov uses builtin from [tawny] -# FUNCTION: DESCRIPTION: -# getCenterRob Extracts the robust estimate for the center -# getCovRob Extracts the robust estimate for the covariance -################################################################################ - - -assetsMeanCov <- -function(x, - method = c("cov", "mve", "mcd", "MCD", "OGK", "nnve", "shrink", "bagged"), - check = TRUE, force = TRUE, baggedR = 100, sigmamu = scaleTau2, alpha = 1/2, - ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Computes robust mean and covariance from multivariate time series - - # Arguments: - # x - a multivariate time series, a data frame, or any other - # rectangular object of assets which can be converted into - # a matrix by the function 'as.matrix'. Optional Dates are - # rownames, instrument names are column names. - # method - Which method should be used to compute the covarinace? - # method = "cov" sample covariance computation - # method = "mve" uses "mve" from [MASS] - # method = "mcd" uses "mcd" from [MASS] - # method = "MCD" uses "MCD" from [robustbase] - # method = "OGK" uses "OGK" from [robustbase] - # method = "nnve" uses "nnve" from [covRobust] - # method = "shrink" uses "shrinkage" from [corpcor] - # method = "bagged" uses "bagging" [corpcor] - # alpha - MCD: numeric parameter controlling the size of the subsets - # over which the determinant is minimized, i.e., alpha*n observations - # are used for computing the determinant. Allowed values are between - # 0.5 and 1 and the default is 0.5. - # sigma.mu - OGK: a function that computes univariate robust location - # and scale estimates. By default it should return a single numeric - # value containing the robust scale (standard deviation) estimate. - # When mu.too is true, sigmamu() should return a numeric vector of - # length 2 containing robust location and scale estimates. See - # scaleTau2, s_Qn, s_Sn, s_mad or s_IQR for examples to be used as - # sigmamu argument. - - # Note: - # The output of this function can be used for portfolio - # optimization. - - # Example: - # DJ = 100 * returns(as.timeSeries(data(DowJones30))) - # DJ = DJ[, c("CAT", "IBM", "GE", "JPM")] - # Sample Covariance: - # assetsMeanCov(DJ, "cov") - # MASS: - # assetsMeanCov(DJ, "mve") - # assetsMeanCov(DJ, "mcd") - # require(robustbase) - # assetsMeanCov(DJ, "MCD") - # assetsMeanCov(DJ, "OGK") - # require(covRobust) - # assetsMeanCov(DJ, "nnve") - - # FUNCTION: - - # Transform Input: - x.mat = as.matrix(x) - - # Do not use: method = match.arg(method) - method = method[1] - N = ncol(x) - assetNames = colnames(x.mat) - - # Attribute Control List: - control = c(method = method[1]) - user = TRUE - - # Compute Classical Covariance: - if (method == "cov") { - # Classical Covariance Estimation: - ans = list(center = colMeans(x.mat), cov = cov(x.mat)) - user = FALSE - } - - # From R Package "robustbase": - if (method == "MCD" | method == "Mcd") { - ans = robustbase::covMcd(x.mat, alpha = alpha, ...) - mu = ans$center - Sigma = ans$cov - user = FALSE - } - if (method == "OGK" | method == "Ogk") { - ans = robustbase::covOGK(x.mat, sigmamu = scaleTau2, ...) - user = FALSE - } - - # [MASS] mve and mcd Routines: - if (method == "mve") { - # require(MASS) - ans = MASS::cov.rob(x = x.mat, method = "mve") - user = FALSE - } - if (method == "mcd") { - # require(MASS) - ans = MASS::cov.rob(x = x.mat, method = "mcd") - user = FALSE - } - - # [corpcor] Shrinkage and Bagging Routines - if (method == "shrink") { - fit = .cov.shrink(x = x.mat, ...) - ans = list(center = colMeans(x.mat), cov = fit) - user = FALSE - } - if (method == "bagged") { - fit = .cov.bagged(x = x.mat, R = baggedR, ...) - ans = list(center = colMeans(x.mat), cov = fit) - control = c(control, R = as.character(baggedR)) - user = FALSE - } - - # Nearest Neighbour Variance Estimation: - if (method == "nnve") { - fit = .cov.nnve(datamat = x.mat, ...) - ans = list(center = colMeans(x.mat), cov = fit$cov) - user = FALSE - } - - # User specified estimator: - if(user) { - fun = match.fun(method[1]) - ans = fun(x.mat, ...) - } - - # Result: - mu = center = ans$center - Sigma = cov = ans$cov - - # Add Size to Control List: - control = c(control, size = as.character(N)) - - # Add Names for Covariance Matrix to Control List: - names(mu) = assetNames - colnames(Sigma) = rownames(Sigma) = colNames = assetNames - - # Check Positive Definiteness: - if (check) { - result = isPositiveDefinite(Sigma) - if(result) { - control = c(control, posdef = "TRUE") - } else { - control = c(control, posdef = "FALSE") - } - } - - # Check Positive Definiteness: - control = c(control, forced = "FALSE") - if (force) { - control = c(control, forced = "TRUE") - if (!result) Sigma = makePositiveDefinite(Sigma) - } - - # Result: - ans = list(center = mu, cov = Sigma, mu = mu, Sigma = Sigma) - attr(ans, "control") = control - - # Return Value: - ans -} - - -################################################################################ - - -.covMeanCov <- -function(x, ...) -{ - # Description: - # Uses sample covariance estimation - - # Arguments: - # x - an object of class timeSeries - - # FUNCTION: - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - - ans = list(center = colMeans(x.mat), cov = cov(x.mat)) - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.mveMeanCov <- -function(x, ...) -{ - # Description: - - # Arguments: - # x - an object of class timeSeries - - # FUNCTION: - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - - ans = MASS::cov.rob(x = x.mat, method = "mve") - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.mcdMeanCov <- -function(x, ...) -{ - # Description: - - # Arguments: - # x - an object of class timeSeries - - # FUNCTION: - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - - ans = MASS::cov.rob(x = x.mat, method = "mcd") - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.studentMeanCov <- -function(x, ...) -{ - # Description: - - # Arguments: - # x - an object of class timeSeries - - # FUNCTION: - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - - ans = MASS::cov.trob(x, ...) - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.MCDMeanCov <- -function(x, alpha = 1/2, ...) -{ - # Description: - - # Arguments: - # x - an object of class timeSeries - - # FUNCTION: - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - - ans = robustbase::covMcd(x.mat, alpha = alpha, ...) - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.OGKMeanCov <- -function(x, sigmamu = scaleTau2, ...) -{ - # Description: - - # Arguments: - # x - an object of class timeSeries - - # FUNCTION: - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - - ans = robustbase::covOGK(x.mat, sigmamu = scaleTau2, ...) - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.nnveMeanCov <- -function(x, ...) -{ - # Description: - - # Arguments: - # x - an object of class timeSeries - - # FUNCTION: - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - fit = .cov.nnve(datamat = x.mat, ...) - - ans = list(center = colMeans(x.mat), cov = fit$cov) - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.shrinkMeanCov <- -function(x, ...) -{ - # Description: - - # Arguments: - # x - an object of class timeSeries - - # Note: - # Based on a function borrowed from package corpcor - - # FUNCTION: - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - fit = .cov.shrink(x = x.mat, ...) - - ans = list(center = colMeans(x.mat), cov = fit) - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.baggedMeanCov <- -function(x, baggedR = 100, ...) -{ - # Description: - - # Arguments: - # x - an object of class timeSeries - - # Note: - # Based on a function borrowed from package corpcor - - # FUNCTION: - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - fit = .cov.bagged(x = x.mat, R = baggedR, ...) - - ans = list(center = colMeans(x.mat), cov = fit) - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.arwMeanCov <- -function(x, ...) -{ - # Description: - # Adaptive reweighted estimator for multivariate location and scatter - # with hard-rejection weights and delta = chi2inv(1-d,p) - - # Arguments: - # x - an object of class timeSeries - - # Note: - # Based on a function borrowed from package mvoutlier - - # FUNCTION: - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - fit = .cov.arw(x = x.mat, center = colMeans(x.mat), cov = cov(x),, ...) - - ans = list(center = fit$center, cov = fit$cov) - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.donostahMeanCov <- -function(x, ...) -{ - # Description: - - # Arguments: - # x - an object of class timeSeries - - # Note: - # Based on a function borrowed from package robust - - # Settings: - x.mat = as.matrix(x) - N = ncol(x) - assetNames = colnames(x) - - ans = .cov.donostah(x = x.mat) - names(ans$center) = assetNames - rownames(ans$cov) = colnames(ans$cov) = assetNames - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.bayesSteinMeanCov <- -function(x, ...) -{ - # Description: - - # Arguments: - # x - an object of class timeSeries - - # Note: - # Based on a function written by Alexios Ghalanos - - # Bayes Stein estimator - # Alexios Ghalanos 2008 - # alexios at 4dscape.com - # This function encapsulates an example of shrinking the returns - # and covariance using Bayes-Stein shrinkage as described in - # Jorion, 1986. - - # Settings: - data <- getDataPart(x) - mu <- as.matrix(apply(data,2, FUN = function(x) mean(x))) - S <- cov(data) - k <- dim(data)[2] - n <- dim(data)[1] - one <- as.matrix(rep(1, k)) - a <- solve(S, one) - - # Constant non informative prior - mu.prior <- one * as.numeric(t(mu) %*% a/t(one) %*% a) - S.inv <- solve(S) - d <- t(mu-mu.prior) %*% S.inv %*% (mu-mu.prior) - d <- as.numeric(d) - lambda <- (k+2) / d - w <- lambda / (n+lambda) - mu.pred <- (1-w) * mu + w * mu.prior - - wc1 <- 1 / (n+lambda) - wc2 <- lambda*(n-1) / (n*(n+lambda)*(n-k-2)) - wc2 <- wc2 / as.numeric(t(one) %*% a) - V.post <- wc1 * S + wc2 * one %*% t(one) - V.pred <- S + V.post - sigma.post <- sqrt(diag(V.post)) - sigma.pred <- sqrt(diag(V.pred)) - - result <- list( - mu = mu, mu.prior = mu.prior, mu.predict = mu.pred, - V = S, V.post = V.post, V.pred = V.pred, Sigma = sqrt(diag(S)), - Sigma.post = sigma.post, Sigma.predict = sigma.pred) - - ans = list(center = result$mu.pred[,1], cov = result$V.pred) - names(ans$center) = colnames(x) - rownames(ans$cov) = colnames(ans$cov) = colnames(x) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.ledoitWolfMeanCov <- -function(x, ...) -{ - # Description: - # Perform shrinkage on a sample covariance towards a biased covariance - - # Arguments: - # x - an object of class timeSeries - - # Details: - # This performs a covariance shrinkage estimation as specified in - # Ledoit and Wolf. Using within the larger framework only requires - # using the getCorFilter.Shrinkage function, which handles the work - # of constructing a shrinkage estimate of the covariance matrix of - # returns (and consequently its corresponding correlation matrix). - - # Note: - # Based on a function borrowed from package tawny - # Author: Brian Lee Yung Rowe - - # Settings: - data = getDataPart(x) - center = colMeans(data) - cov = .cov.shrink.tawny(data, ...) - - ans = list(center = center, cov = cov) - names(ans$center) = colnames(x) - rownames(ans$cov) = colnames(ans$cov) = colnames(x) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -.rmtMeanCov <- -function(x, ...) -{ - # Description: - # Perform Random Matrix Theory on correlation matrix - - # Arguments: - # x - an object of class timeSeries - - # Author: - # tawnyBrian Lee Yung Rowe - - # Note: - # Based on a function borrowed from package tawny - # Author: Brian Lee Yung Rowe - - # FUNCTION: - - # Settings: - data = getDataPart(x) - center = colMeans(data) - cor = .filter.RMT(data, trace = FALSE, doplot = FALSE) - - g = colSds(data) - N = length(g) - cov = 0*cor - for (i in 1:N) - for (j in i:N) - cov[i,j] = cov[j,i] = g[i]*cor[i,j]*g[j] - - ans = list(center = center, cov = cov) - names(ans$center) = colnames(x) - rownames(ans$cov) = colnames(ans$cov) = colnames(x) - - # Return Value: - ans -} - - -################################################################################ - - -getCenterRob <- -function(object) -{ - # Return Value: - object$center -} - - -# ------------------------------------------------------------------------------ - - -getCovRob <- -function(object) -{ - # Return Value: - object$cov -} - - -################################################################################ - diff -Nru fassets-2100.78/R/assetsOutliers.R fassets-2110.79/R/assetsOutliers.R --- fassets-2100.78/R/assetsOutliers.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/assetsOutliers.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,96 +0,0 @@ - -# This library is free software, you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation, either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY, without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR Description. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library, if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - - -################################################################################ -# FUNCTION: DESCRIPTION: -# assetsOutliers Detects outliers in multivariate assets sets -################################################################################ - - -assetsOutliers <- - function (x, center, cov, ...) -{ - # An adapted copy from contributed R package mvoutlier - - # Description: - # Detects outliers in a multivariate set of assets - - # Arguments: - - # Source: - # The code concerned with the outliers is from R package "mvoutliers" - # Moritz Gschwandtner - # Peter Filzmoser - - # References: - # P. Filzmoser, R.G. Garrett, and C. Reimann (2005). - # Multivariate Outlier Detection in Exploration Geochemistry. - # Computers & Geosciences. - - # FUNCTION: - - # Check timeSeries Input: - stopifnot(is.timeSeries(x)) - tS = x - x = as.matrix(x) - - # Critical Values: - n = nrow(x) - p = ncol(x) - if (p <= 10) pcrit = (0.240 - 0.0030 * p)/sqrt(n) - if (p > 10) pcrit = (0.252 - 0.0018 * p)/sqrt(n) - delta = qchisq(0.975, p) - - # Compute Mahalanobis Squared Distances: - d2 = mahalanobis(x, center, cov) - - # Detect Outliers: - d2ord = sort(d2) - dif = pchisq(d2ord, p) - (0.5:n)/n - i = (d2ord >= delta) & (dif > 0) - if (sum(i) == 0) alfan = 0 else alfan = max(dif[i]) - if (alfan < pcrit) alfan = 0 - if (alfan > 0) cn = max(d2ord[n-ceiling(n*alfan)], delta) else cn = Inf - w = d2 < cn - m = apply(x[w, ], 2, mean) - c1 = as.matrix(x - rep(1, n) %*% t(m)) - c = (t(c1) %*% diag(w) %*% c1)/sum(w) - - # Identify Outliers: - outliers = (1:dim(x)[1])[!w] - if (length(outliers) == 0) { - outliers = NA - } else { - names(outliers) = rownames(x)[outliers] - } - - # Compose Result: - ans = list( - center = m, - cov = c, - cor = cov2cor(c), - quantile = cn, - outliers = outliers, - series = tS[outliers, ]) - - # Return Value: - ans -} - - -################################################################################ - diff -Nru fassets-2100.78/R/assetsPfolio.R fassets-2110.79/R/assetsPfolio.R --- fassets-2100.78/R/assetsPfolio.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/assetsPfolio.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,413 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR Description. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - Diethelm Wuertz, GPL -# 2007 - Rmetrics Foundation, GPL -# Diethelm Wuertz -# for code accessed (or partly included) from other sources: -# see Rmetric's copyright and license files - - -################################################################################ -# FUNCTION: DESCRIPTION: -# pfolioVaR Computes VaR for a portfolio of assets -# pfolioCVaR Computes CVaR for a portfoluio of assets -# pfolioCVaRplus Computes CVaR-Plus for a portfolio of assets -# lambdaCVaR Computes CVaR's atomic split value lambda -# FUNCTION: DESCRIPTION: -# pfolioMaxLoss Computes maximum loss for a portfolio -# pfolioReturn Computes return series for a portfolio -# pfolioTargetReturn Computes target return for a portfolio -# pfolioTargetRisk Computes target risk for a portfolio -# pfolioHist Plots a histogram of portfolio returns -################################################################################ - - -pfolioVaR <- - function(x, weights = NULL, alpha = 0.05) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Compute Value-at-Risk for a portfolio of assets - - # Arguments: - # x - a time series, data.frame or any other rectangular object - # of assets which can be written as a matrix object - # weights - a numeric vector of weights - # alpha - the confidence level - - # FUNCTION: - - # Transform: - x = as.matrix(x) - - # Compute Portfolio VaR: - if (is.null(weights)) - weights = rep(1/dim(x)[[2]], dim(x)[[2]]) - n = dim(x)[1] - x = apply(t(t(x) * weights), 1, sum) - n.alpha = max(floor(n * alpha)) - ans = as.vector(sort(x)[n.alpha]) - names(ans) = "VaR" - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -pfolioCVaRplus <- - function(x, weights = NULL, alpha = 0.05) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Compute Value-at-Risk Plus for a portfolio of assets - - # Arguments: - # x - a time series, data.frame or any other rectangular object - # of assets which can be written as a matrix object - # weights - a numeric vector of weights - # alpha - the confidence level - - # FUNCTION: - - # Transform: - x = as.matrix(x) - - # Compute Portfolio CVaRplus: - if (is.null(weights)) { - weights = rep(1/dim(x)[[2]], dim(x)[[2]]) - } - n = dim(x)[1] - x = apply(t(t(x) * weights), 1, sum) - n.alpha = max(1, floor(n * alpha)-1) - ans = as.vector(mean(sort(x)[1:n.alpha])) - names(ans) = "CVaRplus" - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -pfolioCVaR <- - function(x, weights = NULL, alpha = 0.05) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Compute Conditional Value-at-risk for a portfolio of assets - - # Arguments: - # x - a time series, data.frame or any other rectangular object - # of assets which can be written as a matrix object - # weights - a numeric vector of weights - # alpha - the confidence level - # lambda - split value - - # FUNCTION: - - # Transform: - data = as.matrix(x) - - # Input Data: - if (is.null(weights)) { - weights = rep(1/dim(data)[[2]], dim(data)[[2]]) - } - n = dim(data)[1] - Rp = apply(t(t(data)*weights), 1, sum) - - # Sort the Portfolio returns Y - sorted = sort(Rp) - - # Compute Portfolio VaR: - n.alpha = floor(n*alpha) - VaR = sorted[n.alpha] - - # Compute Portfolio CVaRplus: - n.alpha = max(1, floor(n*alpha)-1) - CVaRplus = mean(sorted[1:n.alpha]) - - # Compute Portfolio CVaR: - lambda = 1 - floor(n*alpha)/(n*alpha) - ans = as.vector(lambda*VaR + (1-lambda)*CVaRplus) - names(ans) = "CVaR" - attr(ans, "control") = c(CVaRplus = CVaRplus, lambda = lambda) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -lambdaCVaR <- - function(n, alpha = 0.05) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Computes CVaR's atomic split value lambda - - # Arguments: - # n - the number of oberservations - # alpha - the confidence interval - - # FUNCTION: - - # Compute CVaR lambda: - lambda = 1 - floor(alpha * n) / (alpha * n) - names(lambda) = "lambda" - - # Return Value: - lambda -} - - -################################################################################ - - -pfolioMaxLoss <- - function(x, weights = NULL) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Computes maximum loss for a portfolio of assets - - # Arguments: - # x - a timeSeries, data.frame or any other rectangular object - # of assets which can be written as a matrix object - # weights - the vector of weights - # alpha - the confidence level - - # FUNCTION: - - # Transform: - x = as.matrix(x) - - # Compute MaxLoss [MinReturn]: - if (is.null(weights)) { - weights = rep(1/dim(x)[[2]], dim(x)[[2]]) - } - x = apply(t(t(x)*weights), 1, sum) - ans = min(x) - - # Return Value: - ans -} - - -# ------------------------------------------------------------------------------ - - -pfolioReturn <- - function(x, weights = NULL) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Computes return value of a portfolio - - # Arguments: - # x - a timeSeries, data.frame or any other rectangular object - # of assets which can be written as a matrix object - # weights - the vector of weights - - # FUNCTION: - - # Time Series Object ? - if ( class(x) == "timeSeries" ) { - TS = TRUE - data = x[, 1] - } else { - TS = FALSE - } - - # Compute Portfolio Returns: - x = as.matrix(x) - if (is.null(weights)) - weights = rep(1/dim(x)[[2]], dim(x)[[2]]) - ans = t(t(apply(t(t(x)*weights), 1, sum))) - colnames(ans) = "Portfolio" - - # Time Series Object: - if (TS) { - series(data) = ans - } else { - data = ans - } - - # Return Value: - data -} - - -# ------------------------------------------------------------------------------ - - -pfolioTargetReturn <- - function(x, weights = NULL) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Computes return value of a portfolio - - # Arguments: - # x - a timeSeries, data.frame or any other rectangular object - # of assets which can be written as a matrix object - # weights - the vector of weights - - # FUNCTION: - - # Transform: - x = as.matrix(x) - - # Compute Portfolio Returns: - ans = mean(pfolioReturn(x = x, weights = weights)) - - # Return Value: - names(ans) = "TargetReturn" - ans -} - - -# ------------------------------------------------------------------------------ - - -pfolioTargetRisk <- - function(x, weights = NULL) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Computes risk from covariance matrix of a portfolio - - # Arguments: - # x - a timeSeries, data.frame or any other rectangular object - # of assets which can be written as a matrix object - # weights - the vector of weights - - # FUNCTION: - - # Transform: - x = as.matrix(x) - - # Compute Portfolio Returns: - if (is.null(weights)) - weights = rep(1/dim(x)[[2]], dim(x)[[2]]) - ans = as.vector(sqrt(weights %*% cov(x) %*% weights)) - - # Return Value: - names(ans) = "TargetRisk" - ans -} - - - -# ------------------------------------------------------------------------------ - - -pfolioHist <- - function(x, weights = NULL, alpha = 0.05, range = NULL, details = TRUE, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Plots a histogram of the returns of a portfolio - - # Arguments: - # x - a timeSeries, data.frame or any other rectangular object - # of assets which can be written as a matrix object - # weights - the vector of weights - - # FUNCTION: - - # Transform: - x = as.matrix(x) - - # Suppress Warnings: - opt = options() - options(warn = -1) - - # Plot Portfolio Returns: - Returns = pfolioReturn(x = x, weights = weights) - if (is.null(range)) { - lim = 1.05 * pfolioMaxLoss(x = x, weights = weights)[[1]] - xlim = c(lim, -lim) - } else { - xlim = range - } - Histogram = hist(Returns, xlim = xlim, xlab = "Portfolio Return %", - probability = TRUE, col = "steelblue4", border = "white", ...) - r = seq(xlim[1], xlim[2], length = 201) - lines(r, dnorm(r, mean = mean(Returns), sd = sd(Returns)), ...) - - points(Returns, rep(0, length(Returns)), pch = 20, - col = "orange", cex = 1.25) - - # Add VaR, CVaRplus and MaxLoss: - V1 = pfolioVaR(x = x, weights = weights, alpha = alpha)[[1]] - abline(v = V1, col = "blue", ...) - V2 = pfolioCVaRplus(x = x, weights = weights, alpha = alpha)[[1]] - abline(v = V2, col = "red", ...) - V3 = pfolioMaxLoss(x = x, weights = weights)[[1]] - abline(v = V3, col = "green", ...) - V4 = as.vector(mean(Returns))[[1]] - V5 = as.vector(sd(Returns))[[1]] - - yt = max(density(Returns)$y) - - text(V1, yt, as.character(round(V1, 2)), cex = 0.75, col = "orange") - text(V2, yt, as.character(round(V2, 2)), cex = 0.75, col = "orange") - text(V3, yt, as.character(round(V3, 2)), cex = 0.75, col = "orange") - text(V4, yt, as.character(round(V4, 2)), cex = 0.75, col = "orange") - - yt = 0.95 * yt - text(V1, yt, "VaR", cex = 0.75, col = "orange") - text(V2, yt, "CVaR+", cex = 0.75, col = "orange") - text(V3, yt, "maxLoss", cex = 0.75, col = "orange") - text(V4, yt, "Mean", cex = 0.75, col = "orange") - - # Result: - options(opt) - ans = list(VaR = V1, VaRplus = V2, maxLoss = V3, mean = V4, sd = V5) - if (details) { - cat("\nVaR: ", V1) - cat("\nVaRplus: ", V2) - cat("\nmax Loss: ", V3) - cat("\nMean: ", V4) - cat("\nStDev: ", V5) - cat("\n") - } - - # Return Value: - invisible(ans) -} - - -################################################################################ - diff -Nru fassets-2100.78/R/assetsSelect.R fassets-2110.79/R/assetsSelect.R --- fassets-2100.78/R/assetsSelect.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/assetsSelect.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,121 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - - -################################################################################ -# FUNCTION: DESCRIPTION: -# assetsSelect Selects similar or dissimilar assets -# .hclustSelect Selects due to hierarchical clustering -# .kmeansSelect Selects due to k-means clustering -################################################################################ - - -assetsSelect = - function(x, method = c("hclust", "kmeans"), control = NULL, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Clusters a set of assets - - # Arguments: - # method - which algorithm should be used? - # hclust - Hierarchical clustering on a set of dissimilarities - # kmeans - k-means clustering on a data matrix - - # FUNCTION: - - # Selection: - # do not method = match.arg(method) to allow for user specified clustering - method = method[1] - - # Transform to matrix: - if (class(x) == "timeSeries") { - x = as.matrix(x) - } - - # Compose Function: - fun = paste(".", method, "Select", sep = "") - FUN = match.fun(fun) - - # Cluster: - ans = FUN(x, control, ...) - - # Return Value: - ans -} - - -################################################################################ - - -.hclustSelect <- - function(x, control = NULL, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Hierarchical Clustering - - # FUNCTION: - - # Method: - if (is.null(control)) - control = c(measure = "euclidean", method = "complete") - measure = control[1] - method = control[2] - - # hclust: - ans = hclust(dist(t(x), method = measure), method = method, ...) - class(ans) = c("list", "hclust") - - # Return Value: - ans -} - - -################################################################################ - - -.kmeansSelect <- - function(x, control = NULL, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # kmeans Clustering - - # Note: - # centers must be specified by the user! - - # FUNCTION: - - # Method: - if (is.null(control)) - control = c(centers = 5, algorithm = "Hartigan-Wong") - centers = as.integer(control[1]) - algorithm = control[2] - - # kmeans: - ans = kmeans(x = t(x), centers = centers, algorithm = algorithm, ...) - class(ans) = c("list", "kmeans") - - # Return Value: - ans -} - - -################################################################################ \ No newline at end of file diff -Nru fassets-2100.78/R/assetsSim.R fassets-2110.79/R/assetsSim.R --- fassets-2100.78/R/assetsSim.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/assetsSim.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,89 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - - -################################################################################ -# FUNCTION: DESCRIPTION: -# assetsSim Simulates a set of artificial assets -################################################################################ - - -assetsSim <- - function(n, dim = 2, model = list(mu = rep(0, dim), Omega = diag(dim), - alpha = rep(0, dim), df = Inf), assetNames = NULL) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Simulates a multivariate set of asset log-returns distributed - # according to a Normal, skew-Normal, or skew Student-t Distribution - # Function. - - # Arguments: - # n - the number of data records to be simulated - # dim - the dimension number, i.e. the number of assets to be simulated - # model - a list with the model parameters: - # mu - the numeric vector of mean values of each asset time series - # Omega - the covariance matrix of assets - # alpha - the skewness vector - # df - the degrees of freedom, a measures for the kurtosis - # assetNames - a string vector of asset names, by default NULL - # which creates asset names as "V1", "V2", ..., "Vd", where - # d denotes the dimension - - # Notes: - # Requires function "msn.quantities" from R's GPL licensed - # contributed package "sn", (C) 1998-2004 A. Azzalini. - # The model can also be the value returned by model slot from - # function assetsFit(). - - # Example: - # assetsSim(n=25) - # assetsSim(n=25, assetNames = c("RETURN-1", "RETURN-2") - # assetsSim(n=25, list(mu=c(0,0), Omega=diag(2), alpha=c(0,0), df=4)) - - # FUNCTION: - - # Dimensions: - d = length(model$alpha) - if ( length(model$mu) != d | any(dim(model$Omega) != c(d, d))) - stop("dimensions of arguments do not match") - - # Adapted from contributed R package "sn:rmsn" - Z = msn.quantities(model$mu, model$Omega, model$alpha) - y = matrix(rnorm(n * d), n, d) %*% chol(Z$Psi) - abs.y0 = matrix(rep(abs(rnorm(n)), d), ncol = d) - z = Z$delta * t(abs.y0) + sqrt(1 - Z$delta^2) * t(y) - - # Select: - if (model$df == Inf) { - ans = t(model$mu + Z$omega * z) - } else { - x = rchisq(n, model$df)/model$df - z = t(model$mu + Z$omega * z) - ans = t(model$mu + t(sqrt(x) * z)) - } - - # Dimnames: - dimnames(ans)[[2]] = assetNames - - # Return Value: - as.data.frame(ans) -} - - -################################################################################ - diff -Nru fassets-2100.78/R/assetsTest.R fassets-2110.79/R/assetsTest.R --- fassets-2100.78/R/assetsTest.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/assetsTest.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,191 +0,0 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - - -################################################################################ -# FUNCTION: DESCRIPTION: -# assetsTest Tests for multivariate Normal Assets -# .mvshapiroTest Multivariate Shapiro Test -# .mvenergyTest Multivariate E-Statistic (Energy) Test -################################################################################ - - -assetsTest = -function(x, method = c("shapiro", "energy"), Replicates = 100, -title = NULL, description = NULL) -{ - # Description: - # Tests for multivariate Normal Assets - - # Example: - # .mvnormTest(x = assetsSim(100)) - # .mvnormTest(x = assetsSim(100), method = "e", Replicates = 99) - - # FUNCTION: - - # Test: - method = match.arg(method) - if (method == "shapiro") { - test = .mvshapiroTest(x) - } - if (method == "energy") { - test = .mvenergyTest(x, Replicates = Replicates) - } - - # Return Value: - test -} - - -# ------------------------------------------------------------------------------ - - -.mvshapiroTest = -function(x, title = NULL, description = NULL) -{ - # Description: - # Computes Shapiro's normality test for multivariate variables - - # Note: - # Reimplemented function, doesn't require the contributed R package - # mvnormtest - - # Author: - # Slawomir Jarek - # License: GPL - - # Source: - # Package: mvnormtest - # Version: 0.1-6 - # Date: 2005-04-02 - # Title: Normality test for multivariate variables - # Author: Slawomir Jarek - # Maintainer: Slawomir Jarek - # Description: Generalization of shapiro-wilk test for - # multivariate variables. - - # Example: - # .mvshapiroTest(x = assetsSim(100)) - - # FUNCTION: - - # Transform: - U = t(as.matrix(x)) - - # Test: - n = ncol(U) - if (n < 3 || n > 5000) stop("sample size must be between 3 and 5000") - rng = range(U) - rng = rng[2] - rng[1] - if (rng == 0) - stop("all `U[]' are identical") - Us = apply(U, 1, mean) - R = U-Us - M.1 = solve(R %*% t(R), tol = 1e-18) - Rmax = diag(t(R) %*% M.1 %*% R) - C = M.1 %*% R[, which.max(Rmax)] - Z = t(C) %*% U - test = shapiro.test(Z) - names(test$p.value) = "" - class(test) = "list" - - # Add title and description: - if (is.null(title)) title = "Multivariate Shapiro Test" - if (is.null(description)) description = description() - - # Return Value: - new("fHTEST", - call = match.call(), - data = list(x = x), - test = test, - title = title, - description = description) -} - - -# ------------------------------------------------------------------------------ - - -.mvenergyTest = -function(x, Replicates = 99, title = NULL, description = NULL) -{ - # Description: - # Computes E-statistics test for multivariate variables - - # Note: - # Reimplemented function, doesn't require the contributed - # R package energy, we only use the C Program here. - - # Source: - # Maria L. Rizzo and - # Gabor J. Szekely - # License: GPL 2.0 or later - - # Example: - # .mvenergyTest(x = assetsSim(100), 99) - - # FUNCTION: - - # Transform: - if (class(x) == "timeSeries") x = series(x) - x = as.matrix(x) - - # Test: - R = Replicates - - # RVs: - n <- nrow(x) - d <- ncol(x) - ran.gen = function(x, y) return(matrix(rnorm(n * d), nrow = n, ncol = d)) - - # Parametric Mini Boot: - strata = rep(1, n) - n <- nrow(x) - temp.str <- strata - strata <- tapply(1:n, as.numeric(strata)) - t0 <- .mvnorm.e(x) - lt0 <- length(t0) - t.star <- matrix(NA, sum(R), lt0) - pred.i <- NULL - for(r in 1:R) t.star[r, ] <- .mvnorm.e(ran.gen(x, NULL)) - - # Result: - test <- list( - statistic = c("E-Statistic" = t0), - p.value = 1 - mean(t.star < t0), - method = "Energy Test", - data.name = paste("x, obs ", n, ", dim ", d, ", reps ", R, sep = "")) - names(test$p.value) = "" - class(test) = "list" - - # Add: - if (is.null(title)) title = test$method - if (is.null(description)) description = description() - - # Return Value: - new("fHTEST", - call = match.call(), - data = list(x = x), - test = test, - title = title, - description = description) -} - - -################################################################################ - - - diff -Nru fassets-2100.78/R/builtin-ecodist.R fassets-2110.79/R/builtin-ecodist.R --- fassets-2100.78/R/builtin-ecodist.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/builtin-ecodist.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,276 @@ + +# This library is free software, you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation, either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY, without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR Description. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library, if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: +# distance +################################################################################ + + +# Rmetrics: +# Note that ecodist is not available on Debian as of 2009-09-28. +# To run these functions under Debian/Rmetrics we have them +# implemented here as a builtin. + + +# Package: ecodist +# Version: 1.2.2 +# Date: 2008-12-15 +# Title: Dissimilarity-based functions for ecological analysis +# Author: Sarah Goslee and Dean Urban +# Maintainer: Sarah Goslee +# Depends: stats +# Description: Dissimilarity-based analysis functions including +# ordination and Mantel test functions, intended for use with +# spatial and community data. +# License: GPL version 2 or newer +# Packaged: Mon Dec 15 09:01:37 2008 + + +# distance <- +.ecodist <- +function(x, method="euclidean") +{ + # calculates similarity and dissimilarity coefficients + # as described in Legendre and Legendre 1998 + # returns lower-triangle + ### + # Sarah Goslee + # 2 March 2006 + # revised 31 March 2008 + # bug-fix 15 December 2008 + ### + # uses clever matrix math to calculate the pieces needed + # by dissimilarity matrices, to make it easy to add new + # indices. + ### + # to add a new metric: + # add it to the commented list below + # add it to the end of the METHODS <- c(...) list + # add the code at the appropriate point at the bottom of + # the function + + # 1: euclidean + # 2: bray-curtis + # 3: manhattan + # 4: mahalanobis + # 5: jaccard + # 6: simple difference + # 7: sorensen + + pairedsum <- function(x) + { + ### paired sums + ### returns an N by N by P matrix containing each + ### combination of + N <- nrow(x) + P <- ncol(x) + A <- numeric(N * N * P) + A <- .C("psum", + as.double(as.vector(t(x))), + as.integer(N), + as.integer(P), + A = as.double(A), + PACKAGE = "fAssets")$A + + A <- array(A, dim=c(N, N, P)) + A + } + + paireddiff <- function(x) + { + ### paired differences + N <- nrow(x) + P <- ncol(x) + A <- numeric(N * N * P) + A <- .C("pdiff", + as.double(as.vector(t(x))), + as.integer(N), + as.integer(P), + A = as.double(A), + PACKAGE = "fAssets")$A + + A <- array(A, dim=c(N, N, P)) + A + } + + jointpresence <- function(x) + { + ### joint count of presences + N <- nrow(x) + P <- ncol(x) + A <- numeric(N * N * P) + A <- .C("jpres", + as.double(as.vector(t(x))), + as.integer(N), + as.integer(P), + A = as.double(A), + PACKAGE = "fAssets")$A + + A <- array(A, dim=c(N, N, P)) + A + } + + jointabsence <- function(x) + { + ### joint count of absences + N <- nrow(x) + P <- ncol(x) + A <- numeric(N * N * P) + A <- .C("jabs", + as.double(as.vector(t(x))), + as.integer(N), + as.integer(P), + A = as.double(A), + PACKAGE = "fAssets")$A + + A <- array(A, dim=c(N, N, P)) + A + } + + firstonly <- function(x) + { + ### present only in first sample + N <- nrow(x) + P <- ncol(x) + A <- numeric(N * N * P) + A <- .C("jfirst", + as.double(as.vector(t(x))), + as.integer(N), + as.integer(P), + A = as.double(A), + PACKAGE = "fAssets")$A + + A <- array(A, dim=c(N, N, P)) + A + } + + secondonly <- function(x) + { + ### present only in second sample + N <- nrow(x) + P <- ncol(x) + A <- numeric(N * N * P) + A <- .C("jsec", + as.double(as.vector(t(x))), + as.integer(N), + as.integer(P), + A = as.double(A), + PACKAGE = "fAssets")$A + + A <- array(A, dim=c(N, N, P)) + A + } + + x <- as.matrix(x) + + ## code borrowed from dist() + METHODS <- c( + "euclidean", "bray-curtis", "manhattan", + "mahalanobis", "jaccard", "difference", + "sorensen") + + method <- pmatch(method, METHODS) + if (is.na(method)) + stop("invalid distance method") + if (method == -1) + stop("ambiguous distance method") + N <- nrow(x) + P <- ncol(x) + + + if(method == 1) + { + # Euclidean distance + A <- paireddiff(x) + D <- sqrt(apply(A, 1:2, function(x)sum(x * x))) + } + + if(method == 2) + { + # Bray-Curtis distance + A <- paireddiff(x) + A <- apply(A, 1:2, function(x)sum(abs(x))) + B <- pairedsum(x) + B <- apply(B, 1:2, sum) + D <- A / B + } + + if(method == 3) + { + # unstandardized manhattan distance + A <- paireddiff(x) + D <- apply(A, 1:2, function(x)sum(abs(x))) + } + + if(method == 4) + { + # pairwise Mahalanobis distance + # same as mahal() + icov <- solve(cov(x)) + A <- paireddiff(x) + A1 <- apply(A, 1, function(z)(z %*% icov %*% t(z))) + D <- A1[seq(1, N*N, by=(N+1)), ] + } + + if(method == 5) + { + # Jaccard distance + A <- jointpresence(x) + A <- apply(A, 1:2, sum) + B <- firstonly(x) + B <- apply(B, 1:2, sum) + C <- secondonly(x) + C <- apply(C, 1:2, sum) + D <- 1 - A / (A + B + C) + } + + if(method == 6) + { + # simple difference, NOT symmetric + D <- paireddiff(x)[,,1, drop=TRUE] + } + + if(method == 7) + { + # Sorensen distance + A <- jointpresence(x) + A <- apply(A, 1:2, sum) + B <- firstonly(x) + B <- apply(B, 1:2, sum) + C <- secondonly(x) + C <- apply(C, 1:2, sum) + D <- 1 - (2*A) / (2*A + B + C) + } + + + ## Make the results lower triangular + D <- D[col(D) < row(D)] + + ## give the results attributes similar to dist() + attr(D, "Size") <- N + attr(D, "Labels") <- rownames(x) + attr(D, "Diag") <- FALSE + attr(D, "Upper") <- FALSE + attr(D, "method") <- METHODS[method] + class(D) <- "dist" + D +} + + +################################################################################ + diff -Nru fassets-2100.78/R/plot-pairs.R fassets-2110.79/R/plot-pairs.R --- fassets-2100.78/R/plot-pairs.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/plot-pairs.R 2010-05-05 07:44:43.000000000 +0000 @@ -1,299 +1,300 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - - -################################################################################ -# FUNCTION: DESCRIPTION: -# assetsPairsPlot Displays pairs of scatterplots of assets -# assetsCorgramPlot Displays pairwise correlations between assets -# assetsCorTestPlot Displays and tests pairwise correlations -# assetsCorImagePlot Displays an image plot of a correlations -################################################################################ - - -assetsPairsPlot <- - function(x, labels = TRUE, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Displays pairs of scatterplots of individual assets - - # Arguments: - # x - a timeSeries object or any other rectangular object - # which can be transformed by the function as. matrix - # into a numeric matrix. - # labels - a logical flag. Should default labels be printed? - # Not implemented. - - # Example: - # x = as.timeSeries(data(LPP2005REC))[, 1:6] - # assetsPairsPlot(x) - - # FUNCTION: - - # Settings: - x = as.matrix(x) - - # Pairs Plot: - # Suppress warnings for tick = 0 in ... - warn = options()$warn - options(warn = -1) - pairs(x, ...) - options(warn = warn) - - # Return Value: - invisible() -} - - -# ------------------------------------------------------------------------------ - - -assetsCorgramPlot <- - function(x, labels = TRUE, method = c( "pie", "shade"), ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Displays correlations between assets - - # Arguments: - # x - a timeSeries object or any other rectangular object - # which can be transformed by the function as. matrix - # into a numeric matrix. - # labels - a logical flag. Should default labels be printed? - # Not implemented. - - # Example: - # x = as.timeSeries(data(LPP2005REC))[, 1:6] - # assetsCorgramPlot(x, method = "pie") - # assetsCorgramPlot(x, method = "shade") - # assetsCorgramPlot(x, method = "hist") # ... has a bug, check - - # FUNCTION: - - # Settings: - method <<- match.arg(method) - stopifnot(is.timeSeries(x)) - x = series(x) - - # Internal Function: - .panel.lower = function(x, y, ...) - { - if (method[1] == "pie") { - .panel.pie(x, y, ...) - .panel.pts(x, y, ...) - } else if (method[1] == "shade") { - .panel.shade(x, y, ...) - .panel.pts(x, y, ...) - } else if (method[1] == "hist") { - .panel.shade(x, y, ...) - .panel.hist(x, y, ...) - } - } - .panel.upper = function(x, y, ...) - { - .panel.ellipse(x, y, ...) - } - - # Plot Corellogram - Pies and Ellipses: - .corrgram(x, - lower.panel = .panel.lower, - upper.panel = .panel.upper, - text.panel = .panel.txt, ...) - - # Return Value: - invisible() -} - - -# ------------------------------------------------------------------------------ - - -assetsCorTestPlot <- - function(x, labels = TRUE, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Displays and tests pairwise correlations of assets - - # Arguments: - # x - a timeSeries object or any other rectangular object - # which can be transformed by the function as. matrix - # into a numeric matrix. - # labels - a logical flag. Should default labels be printed? - # Not implemented. - - # Example: - # x = as.timeSeries(data(LPP2005REC))[, 1:6] - # assetsCorTestPlot(x) - - # FUNCTION: - - # Settings: - x = as.matrix(x) - - # Upper Plot Function: - cortestPanel <- - function(x, y, cex, col, ...) - { - if (missing(col)) col = NULL - usr = par("usr"); on.exit(par(usr)) - par(usr = c(0, 1, 0, 1)) - r = abs(cor(x, y)) - txt = format(c(r, 0.123456789), digits = 3)[1] - test = cor.test(x, y) - Signif = symnum(test$p.value, corr = FALSE, na = FALSE, - cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), - symbols = c("*** ", "** ", "* ", ". ", " ")) - text(0.5, 0.5, txt, cex = 1, col = NULL, ...) - text(0.8, 0.8, Signif, cex = 1.5, col = col, ...) - } - - # Lower Plot Function: - lowessPanel = - function (x, y, ...) - { - points(x, y, ...) - ok = is.finite(x) & is.finite(y) - if (any(ok)) lines(lowess(x[ok], y[ok]), col = "brown") - } - - # Plot: - pairs(x, - lower.panel = lowessPanel, - upper.panel = cortestPanel, ...) - - # Return Value: - invisible() -} - - -# ------------------------------------------------------------------------------ - - -assetsCorImagePlot <- - function(x, - labels = TRUE, - show = c("cor", "test"), - use = c("pearson", "kendall", "spearman"), - abbreviate = 3, ...) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Creates an image plot of a correlations - - # Arguments: - # R - data to be evaluated against its own members - - # Details: - # uses relative colors to indicate the strength of the pairwise - # correlation. - - # Authors: - # Sandrine Dudoit, sandrine@stat.berkeley.edu, from "SMA" library - # modified by Peter Carl - # extended by Diethelm Wuertz - - # Example: - # x = as.timeSeries(data(LPP2005REC)) - # assetsCorImagePlot(x[,assetsArrange(x, "hclust")], abbreviate = 5) - - # FUNCTION: - - # Settings: - R = x - - # Match Arguments: - show = match.arg(show) - use = match.arg(use) - - # Handle Missing Values: - R = na.omit(R, ...) - - # Abbreviate Instrument Names: - Names = colnames(R) = substring(colnames(R), 1, abbreviate) - - # Compute Correlation Matrix: - R = as.matrix(R) - n = NCOL(R) - if (show == "cor") { - corr <- cor(R, method = use) - if (show == "test") { - test = corr*NA - for ( i in 1:n) - for (j in 1:n) - test[i,j] = cor.test(R[,i], R[,j], method = use)$p.value - } - } else if (show == "robust") { - stop("robust: Not Yet Implemented") - } else if (show == "shrink") { - stop("robust: Not Yet Implemented") - } - - - ## compute colors for correlation matrix: - corrMatrixcolors <- function (ncolors) - { - k <- round(ncolors/2) - r <- c(rep(0, k), seq(0, 1, length = k)) - g <- c(rev(seq(0, 1, length = k)), rep(0, k)) - b <- rep(0, 2 * k) - res <- (rgb(r,g,b)) - res - } - - ## Plot Image: - ncolors <- 10*length(unique(as.vector(corr))) - image(x = 1:n, y = 1:n, z = corr[, n:1], - col = corrMatrixcolors(ncolors), - axes = FALSE, main = "", xlab = "", ylab = "", ...) - - # Add Text Values: - if (show == "cor") X = t(corr) else X = t(test) - coord = grid2d(1:n, 1:n) - for (i in 1:(n*n)) { - text(coord$x[i], coord$y[n*n+1-i], - round(X[coord$x[i], coord$y[i]], digits = 2), - col = "white", cex = 0.7) - } - - # Add Axis Labels: - if(labels) { - axis(2, at = n:1, labels = Names, las = 2) - axis(1, at = 1:n, labels = Names, las = 2) - Names = c( - pearson = "Pearson", kendall = "Kendall", spearman = "Spearman") - if (show == "test") Test = "Test" else Test = "" - title( - main = paste(Names[use], " Correlation ", Test, " Image", sep = "")) - mText = paste("Method:", show) - mtext(mText, side = 4, adj = 0, col = "grey", cex = 0.7) - } - - # Add Box: - box() - - # Return Value: - invisible() -} - - -################################################################################ - + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: DESCRIPTION: +# assetsPairsPlot Displays pairs of scatterplots of assets +# assetsCorgramPlot Displays pairwise correlations between assets +# assetsCorTestPlot Displays and tests pairwise correlations +# assetsCorImagePlot Displays an image plot of a correlations +################################################################################ + + +assetsPairsPlot <- + function(x, labels = TRUE, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Displays pairs of scatterplots of individual assets + + # Arguments: + # x - a timeSeries object or any other rectangular object + # which can be transformed by the function as. matrix + # into a numeric matrix. + # labels - a logical flag. Should default labels be printed? + # Not implemented. + + # Example: + # x = as.timeSeries(data(LPP2005REC))[, 1:6] + # assetsPairsPlot(x) + + # FUNCTION: + + # Settings: + x = as.matrix(x) + + # Pairs Plot: + # Suppress warnings for tick = 0 in ... + warn = options()$warn + options(warn = -1) + pairs(x, ...) + options(warn = warn) + + # Return Value: + invisible() +} + + +# ------------------------------------------------------------------------------ + + +assetsCorgramPlot <- + function(x, labels = TRUE, method = c( "pie", "shade"), ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Displays correlations between assets + + # Arguments: + # x - a timeSeries object or any other rectangular object + # which can be transformed by the function as. matrix + # into a numeric matrix. + # labels - a logical flag. Should default labels be printed? + # Not implemented. + # Added in call to .corrgram by DJS, 20/02/2010 + + # Example: + # x = as.timeSeries(data(LPP2005REC))[, 1:6] + # assetsCorgramPlot(x, method = "pie") + # assetsCorgramPlot(x, method = "shade") + # assetsCorgramPlot(x, method = "hist") # ... has a bug, check + + # FUNCTION: + + # Settings: + method <<- match.arg(method) + stopifnot(is.timeSeries(x)) + x = series(x) + + # Internal Function: + .panel.lower = function(x, y, ...) + { + if (method[1] == "pie") { + .panel.pie(x, y, ...) + .panel.pts(x, y, ...) + } else if (method[1] == "shade") { + .panel.shade(x, y, ...) + .panel.pts(x, y, ...) + } else if (method[1] == "hist") { + .panel.shade(x, y, ...) + .panel.hist(x, y, ...) + } + } + .panel.upper = function(x, y, ...) + { + .panel.ellipse(x, y, ...) + } + + # Plot Corellogram - Pies and Ellipses: + .corrgram(x, labels = labels, + lower.panel = .panel.lower, + upper.panel = .panel.upper, + text.panel = .panel.txt, ...) + + # Return Value: + invisible() +} + + +# ------------------------------------------------------------------------------ + + +assetsCorTestPlot <- + function(x, labels = TRUE, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Displays and tests pairwise correlations of assets + + # Arguments: + # x - a timeSeries object or any other rectangular object + # which can be transformed by the function as. matrix + # into a numeric matrix. + # labels - a logical flag. Should default labels be printed? + # Not implemented. + + # Example: + # x = as.timeSeries(data(LPP2005REC))[, 1:6] + # assetsCorTestPlot(x) + + # FUNCTION: + + # Settings: + x = as.matrix(x) + + # Upper Plot Function: + cortestPanel <- + function(x, y, cex, col, ...) + { + if (missing(col)) col = NULL + usr = par("usr"); on.exit(par(usr)) + par(usr = c(0, 1, 0, 1)) + r = abs(cor(x, y)) + txt = format(c(r, 0.123456789), digits = 3)[1] + test = cor.test(x, y) + Signif = symnum(test$p.value, corr = FALSE, na = FALSE, + cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1), + symbols = c("*** ", "** ", "* ", ". ", " ")) + text(0.5, 0.5, txt, cex = 1, col = NULL, ...) + text(0.8, 0.8, Signif, cex = 1.5, col = col, ...) + } + + # Lower Plot Function: + lowessPanel = + function (x, y, ...) + { + points(x, y, ...) + ok = is.finite(x) & is.finite(y) + if (any(ok)) lines(lowess(x[ok], y[ok]), col = "brown") + } + + # Plot: + pairs(x, + lower.panel = lowessPanel, + upper.panel = cortestPanel, ...) + + # Return Value: + invisible() +} + + +# ------------------------------------------------------------------------------ + + +assetsCorImagePlot <- + function(x, + labels = TRUE, + show = c("cor", "test"), + use = c("pearson", "kendall", "spearman"), + abbreviate = 3, ...) +{ + # A function implemented by Diethelm Wuertz + + # Description: + # Creates an image plot of a correlations + + # Arguments: + # R - data to be evaluated against its own members + + # Details: + # uses relative colors to indicate the strength of the pairwise + # correlation. + + # Authors: + # Sandrine Dudoit, sandrine@stat.berkeley.edu, from "SMA" library + # modified by Peter Carl + # extended by Diethelm Wuertz + + # Example: + # x = as.timeSeries(data(LPP2005REC)) + # assetsCorImagePlot(x[,assetsArrange(x, "hclust")], abbreviate = 5) + + # FUNCTION: + + # Settings: + R = x + + # Match Arguments: + show = match.arg(show) + use = match.arg(use) + + # Handle Missing Values: + R = na.omit(R, ...) + + # Abbreviate Instrument Names: + Names = colnames(R) = substring(colnames(R), 1, abbreviate) + + # Compute Correlation Matrix: + R = as.matrix(R) + n = NCOL(R) + if (show == "cor") { + corr <- cor(R, method = use) + if (show == "test") { + test = corr*NA + for ( i in 1:n) + for (j in 1:n) + test[i,j] = cor.test(R[,i], R[,j], method = use)$p.value + } + } else if (show == "robust") { + stop("robust: Not Yet Implemented") + } else if (show == "shrink") { + stop("robust: Not Yet Implemented") + } + + + ## compute colors for correlation matrix: + corrMatrixcolors <- function (ncolors) + { + k <- round(ncolors/2) + r <- c(rep(0, k), seq(0, 1, length = k)) + g <- c(rev(seq(0, 1, length = k)), rep(0, k)) + b <- rep(0, 2 * k) + res <- (rgb(r,g,b)) + res + } + + ## Plot Image: + ncolors <- 10*length(unique(as.vector(corr))) + image(x = 1:n, y = 1:n, z = corr[, n:1], + col = corrMatrixcolors(ncolors), + axes = FALSE, main = "", xlab = "", ylab = "", ...) + + # Add Text Values: + if (show == "cor") X = t(corr) else X = t(test) + coord = grid2d(1:n, 1:n) + for (i in 1:(n*n)) { + text(coord$x[i], coord$y[n*n+1-i], + round(X[coord$x[i], coord$y[i]], digits = 2), + col = "white", cex = 0.7) + } + + # Add Axis Labels: + if(labels) { + axis(2, at = n:1, labels = Names, las = 2) + axis(1, at = 1:n, labels = Names, las = 2) + Names = c( + pearson = "Pearson", kendall = "Kendall", spearman = "Spearman") + if (show == "test") Test = "Test" else Test = "" + title( + main = paste(Names[use], " Correlation ", Test, " Image", sep = "")) + mText = paste("Method:", show) + mtext(mText, side = 4, adj = 0, col = "grey", cex = 0.7) + } + + # Add Box: + box() + + # Return Value: + invisible() +} + + +################################################################################ + diff -Nru fassets-2100.78/R/stats-distance.R fassets-2110.79/R/stats-distance.R --- fassets-2100.78/R/stats-distance.R 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/R/stats-distance.R 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1,314 @@ + +# This library is free software, you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation, either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY, without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR Description. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library, if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +# FUNCTION: +# .corDist +# .kendallDist +# .spearmanDist +# .mutinfoDist +# FUNCTION: +# .euclideanDist +# .maximumDist +# .manhattanDist +# .canberraDist +# .binaryDist +# .minkowskiDist +# FUNCTION: +# .braycurtisDist +# .mahalanobisDist +# .jaccardDist +# .differenceDist +# .sorensenDist +################################################################################ + + +.corDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + x = t(as.matrix(x)) + dist = as.dist(1-cor(x)) + + # Return Value: + dist +} + + +.kendallDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + x = t(as.matrix(x)) + dist = as.dist(1-cor(x, method = "kendall")) + + # Return Value: + dist +} + + +.spearmanDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + x = t(as.matrix(x)) + dist = as.dist(1-cor(x, method = "spearman")) + + # Return Value: + dist +} + + +.mutinfoDist <- +function(x, nbin=10) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # borrowed from R package bioDist and slightly modified + + # Distance: + x <- as.matrix(x) + nc <- ncol(x) + nr <- nrow(x) + clist <- vector("list", length=nr) + for(i in 1:nr) clist[[i]] <- cut(x[i,], breaks=nbin) + ppfun <- function(pp) {pp<-pp[pp>0]; -sum(pp*log(pp ))} + appfun <- function(x,y) { + ppfun(table(x)/nc)+ppfun(table(y)/nc) - ppfun(c(table(x, y)/nc))} + mat = matrix(rep(NA, nr*nr), ncol = nr) + for(i in 1:(nr-1)) { + for(j in (i+1):nr) { + mat[i,j] <- mat[j,i]<- appfun(clist[[i]], clist[[j]]) + } + } + mat = 1 - sqrt(1 - exp(-2*mat)) + colnames(mat) = rownames(mat) = rownames(x) + dist = as.dist(mat) + + # Return Value: + dist +} + + +################################################################################ +# from base R: +# "euclidean", "maximum", "manhattan", +# "canberra", "binary", "minkowski" + + +.euclideanDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + dist = dist(x, "euclidean") + + # Return Value: + dist +} + + +.maximumDist <- +function(x) +{ + # FUNCTION: + + # A function implemented by Diethelm Wuertz + + # Distance: + dist = dist(x, "maximum") + + # Return Value: + dist +} + + +.manhattanDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + dist = dist(x, "manhattan") + + # Return Value: + dist +} + + +.canberraDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + dist = dist(x, "canberra") + + # Return Value: + dist +} + + +.binaryDist <- +function(x) +{ + # Distance: + dist = dist(x, "binary") + + # Return Value: + dist +} + + +.minkowskiDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + dist = dist(x, "minkowski") + + # Return Value: + dist +} + + +################################################################################ + + +# from ecodist: +# "euclidean", "bray-curtis", "manhattan", +# "mahalanobis", "jaccard", "difference" +# "sorensen" + + +.braycurtisDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + dist = .ecodist(x, "bray-curtis") + + # Return Value: + dist +} + + +.mahalanobisDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + dist = .ecodist(x, "mahalanobis") + + # Return Value: + dist +} + + +.jaccardDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + dist = .ecodist(x, "jaccard") + + # Return Value: + dist +} + + +.differenceDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + dist = .ecodist(x, "difference") + + # Return Value: + dist +} + + +.mahalanobisDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + dist = .ecodist(x, "mahalanobis") + + # Return Value: + dist +} + + +.sorensenDist <- +function(x) +{ + # A function implemented by Diethelm Wuertz + + # FUNCTION: + + # Distance: + dist = .ecodist(x, "sorensen") + + # Return Value: + dist +} + + +################################################################################ + + diff -Nru fassets-2100.78/R/zzz.Deprecated.R fassets-2110.79/R/zzz.Deprecated.R --- fassets-2100.78/R/zzz.Deprecated.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/zzz.Deprecated.R 2012-09-24 07:27:22.000000000 +0000 @@ -1,20 +1,42 @@ +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + ################################################################################ # FUNCTION: DESCRIPTION: -# .assetsStats Computes statistics of monthly assets sets +# .assetsStats Computes statistics of monthly assets sets +# FUNCTION: DESCRIPTION: +# .dutchPortfolioData Example Data from Engel's Diploma Thesis +# .usPortfolioData Annual US Economics Portfolio Data +# .sm132PortfolioData Example from Scherer, Martin: Chapter 1.32 +# .worldIndexData A data set of World Indexes +# FUNCTION: SIMILARITY PLOTS: +# fixBinHistogram Returns histogram with fixed bins ################################################################################ -.assetsStats <- +.assetsStats <- function(x) -{ +{ # A function implemented by Diethelm Wuertz # Description: # Computes benchmark statistics for a data set of assets with - # monthly data records. - + # monthly data records. + # Details: # The computed statistics values are: # records - number of records (length of time series) @@ -26,49 +48,49 @@ # maxDD - maximum Drawdown # TUW - Time under Water # mMaxLoss - Monthly maximum Loss - # mVaR - Monthly 99% Value-at-Risk - # mModVaR - Monthly 99% Modified Value-at-Risk + # mVaR - Monthly 99% Value-at-Risk + # mModVaR - Monthly 99% Modified Value-at-Risk # mSharpe - Monthly Sharpe Ratio # mModSharpe - Monthly Modified Sharpe Ratio # skPrice - Skewness/Kurtosis Price # The statistics are implemented based on the formulas from - # "Extreme Metrics". They reflect risk measures as used in + # "Extreme Metrics". They reflect risk measures as used in # the hedge fund software from "www.AlternativeSoft.com". # Arguments: # x - asset data set, a matrix (or vector) where the rows # are numbered by "time", and the columns belong to the # individual assets. Monthly values are expected. - + # Value: # The function returns a data frame with the values of the # 12 statistics for each asset. - + # Reference: # "ExtremeMetrics Software", Help Document, Alternative Software, # March 2003, 4 pages. - + # Example: - + # FUNCTION: - + # If x is a vector, make it a matrix: statistics = 14 if (is.null(dim(x))) { - n = 1 - x = matrix(x, length(x)) + n = 1 + x = matrix(x, length(x)) result = matrix(rep(0, times = statistics), ncol = 1) } else { - n = dim(x)[2] + n = dim(x)[2] result = matrix(rep(0, times = statistics*n), ncol = n) } - - # Give Names to Result Matrix: + + # Give Names to Result Matrix: stat.names = c( "Records", "paMean", "paAve", "paVola", "paSkew", "paKurt", "maxDD", "TUW", "mMaxLoss", "mVaR", "mModVaR", "mSharpe", "mModSharpe", "skPrice") - dimnames(result) = list(stat.names, dimnames(x)[[2]]) + dimnames(result) = list(stat.names, dimnames(x)[[2]]) # Loop over all Assets: for (i in 1:n) { @@ -82,48 +104,43 @@ # Annualized volatility from monthly returns: result[4, i] = annualizedVolatility = sqrt(var(r)) # Annualized skewness from monthly returns: - result[5, i] = annualizedSkewness = skewness(r) + result[5, i] = annualizedSkewness = skewness(r) # Annualized Kurtosis from monthly returns: - result[6, i] = annualizedKurtosis = kurtosis(r) + result[6, i] = annualizedKurtosis = kurtosis(r) # Maximum Drawdown of of monthly returns: result[7, i] = maxDrawdown = max(cummax(cumsum(r)) - cumsum(r)) # Time-Under-Water of monthly returns: - result[8, i] = timeUnderWater = + result[8, i] = timeUnderWater = max(diff(which (diff(cummax(cumsum(r))) != 0))) # Maximum Loss of monthly returns: - result[9, i] = maxMonthlyLoss = min(r) + result[9, i] = maxMonthlyLoss = min(r) # Monthly Value at Risk: zc = 2.33 - result[10, i] = monthlyVaR = annualizedMean - - zc * annualizedVolatility + result[10, i] = monthlyVaR = annualizedMean - + zc * annualizedVolatility # Monthly Modified Value at Risk: - p = 0.99; s = annualizedSkewness; k = annualizedKurtosis + p = 0.99; s = annualizedSkewness; k = annualizedKurtosis zcf = zc + (zc*zc-1)*s/6 + zc*(zc*zc-3)*k/24 + zc*(2*zc*zc-5)*s*s/36 - result[11, i] = monthlyModVaR = annualizedMean - - zcf * annualizedVolatility + result[11, i] = monthlyModVaR = annualizedMean - + zcf * annualizedVolatility # Monthly Sharpe Ratio: - result[12, i] = monthlySharpeRatio = - annualizedMean/annualizedVolatility + result[12, i] = monthlySharpeRatio = + annualizedMean/annualizedVolatility # Monthly Modified Sharpe Ratio: - result[13, i] = monthlyModSharpeRatio = annualizedMean/monthlyModVaR + result[13, i] = monthlyModSharpeRatio = annualizedMean/monthlyModVaR # Skewness Kurtosis Price: - result[14, i] = skewnesskurtosisPrice = annualizedMean * + result[14, i] = skewnesskurtosisPrice = annualizedMean * ( monthlyModVaR/monthlyVaR - 1) } - + # Result: - ans = as.data.frame(round(result, digits = 3)) - + ans = as.data.frame(round(result, digits = 3)) + # Return Value: ans -} +} ################################################################################ -# FUNCTION: DESCRIPTION: -# .dutchPortfolioData Example Data from Engel's Diploma Thesis -# .usPortfolioData Annual US Economics Portfolio Data -# .sm132PortfolioData Example from Scherer, Martin: Chapter 1.32 -# .worldIndexData A data set of World Indexes .dutchPortfolioData = @@ -452,46 +469,40 @@ ################################################################################ -# FUNCTION: SIMILARITY PLOTS: -# fixBinHistogram Returns histogram with fixed bins -################################################################################ - -.hist <- - function (x, nbins) -{ - # A function implemented by Diethelm Wuertz - - # Description: - # Returns histogram with fixed bins - - # FUNCTION: - - # Classes: - nclass = nbins + 1 - n = length(x) - xname = paste(deparse(substitute(x), 500), collapse = "\n") - - # Breaks: - breaks = seq(min(x), max(x), length = nclass) - nB = length(breaks) - h = diff(breaks) - - # Compute Counts: - counts = .C("bincount", as.double(x), as.integer(n), as.double(breaks), - as.integer(nB), counts = integer(nB - 1), right = FALSE, - include = TRUE, naok = FALSE, NAOK = FALSE, DUP = FALSE, - PACKAGE = "base")$counts - dens = counts/(n * h) - mids = 0.5 * (breaks[-1] + breaks[-nB]) - - # Histogram: - r = structure(list(breaks = breaks, counts = counts, intensities = dens, - density = dens, mids = mids, xname = xname, equidist = TRUE), - class = "histogram") -} +## .hist <- +## function (x, nbins) +## { +## # A function implemented by Diethelm Wuertz + +## # Description: +## # Returns histogram with fixed bins + +## # FUNCTION: + +## # Classes: +## nclass = nbins + 1 +## n = length(x) +## xname = paste(deparse(substitute(x), 500), collapse = "\n") + +## # Breaks: +## breaks = seq(min(x), max(x), length = nclass) +## nB = length(breaks) +## h = diff(breaks) + +## # Compute Counts: +## counts = .C("bincount", as.double(x), as.integer(n), as.double(breaks), +## as.integer(nB), counts = integer(nB - 1), right = FALSE, +## include = TRUE, naok = FALSE, NAOK = FALSE, DUP = FALSE, +## PACKAGE = "base")$counts +## dens = counts/(n * h) +## mids = 0.5 * (breaks[-1] + breaks[-nB]) + +## # Histogram: +## r = structure(list(breaks = breaks, counts = counts, intensities = dens, +## density = dens, mids = mids, xname = xname, equidist = TRUE), +## class = "histogram") +## } ################################################################################ - - \ No newline at end of file diff -Nru fassets-2100.78/R/zzz.R fassets-2110.79/R/zzz.R --- fassets-2100.78/R/zzz.R 2009-09-02 20:27:41.000000000 +0000 +++ fassets-2110.79/R/zzz.R 2010-05-05 07:44:43.000000000 +0000 @@ -1,47 +1,44 @@ - -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Library General Public -# License as published by the Free Software Foundation; either -# version 2 of the License, or (at your option) any later version. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Library General Public License for more details. -# -# You should have received a copy of the GNU Library General -# Public License along with this library; if not, write to the -# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, -# MA 02111-1307 USA - -# Copyrights (C) -# for this R-port: -# 1999 - 2008, Diethelm Wuertz, Rmetrics Foundation, GPL -# Diethelm Wuertz -# info@rmetrics.org -# www.rmetrics.org -# for the code accessed (or partly included) from other R-ports: -# see R's copyright and license files -# for the code accessed (or partly included) from contributed R-ports -# and other sources -# see Rmetrics's copyright file - - -################################################################################ - - -.First.lib = -function(lib, pkg) -{ - # Load dll: - library.dynam("fAssets", pkg, lib) -} - - -if(!exists("Sys.setenv", mode = "function")) # pre R-2.5.0, use "old form" - Sys.setenv <- Sys.putenv - - - -################################################################################ - + +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Library General Public +# License as published by the Free Software Foundation; either +# version 2 of the License, or (at your option) any later version. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Library General Public License for more details. +# +# You should have received a copy of the GNU Library General +# Public License along with this library; if not, write to the +# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, +# MA 02111-1307 USA + + +################################################################################ +.First.lib = +function(lib, pkg) +{ + +### # Startup Mesage and Desription: +### MSG <- if(getRversion() >= "2.5") packageStartupMessage else message +### dsc <- packageDescription(pkg) +### if(interactive() || getOption("verbose")) { +### # not in test scripts +### MSG(sprintf("Rmetrics Package %s (%s) loaded.", pkg, dsc$Version)) +### } + + setRmetricsOptions(.x.save = NA) + +} + + +.onLoad <- function(libname, pkgname) setRmetricsOptions(.x.save = NA) + + +if(!exists("Sys.setenv", mode = "function")) # pre R-2.5.0, use "old form" + Sys.setenv <- Sys.putenv + + +################################################################################ + diff -Nru fassets-2100.78/debian/changelog fassets-2110.79/debian/changelog --- fassets-2100.78/debian/changelog 2013-05-05 16:25:27.000000000 +0000 +++ fassets-2110.79/debian/changelog 2013-05-05 16:25:27.000000000 +0000 @@ -1,3 +1,27 @@ +fassets (2110.79-2precise0) precise; urgency=low + + * Compilation for Ubuntu 12.04.2 LTS + + -- Michael Rutter Sun, 05 May 2013 16:16:59 +0000 + +fassets (2110.79-2) unstable; urgency=low + + * debian/control: Set Build-Depends: to current R version + + * (Re-)building with R 3.0.0 + + -- Dirk Eddelbuettel Wed, 03 Apr 2013 05:39:36 -0500 + +fassets (2110.79-1) unstable; urgency=low + + * New upstream release + + * debian/control: Set Build-Depends: to current R version + * debian/control: Change Depends to ${R:Depends} + * debian/control: Set Standards-Version: to current version + + -- Dirk Eddelbuettel Mon, 24 Sep 2012 06:26:04 -0500 + fassets (2100.78-3) unstable; urgency=low * debian/control: Switch (Build-)Depends: from r-cran-vr to r-cran-mass diff -Nru fassets-2100.78/debian/control fassets-2110.79/debian/control --- fassets-2100.78/debian/control 2013-05-05 16:25:27.000000000 +0000 +++ fassets-2110.79/debian/control 2013-05-05 16:25:27.000000000 +0000 @@ -2,13 +2,13 @@ Section: gnu-r Priority: optional Maintainer: Dirk Eddelbuettel -Build-Depends: debhelper (>= 5.0.0), r-base-dev (>= 2.11.0~20100409), cdbs, r-cran-sn, r-cran-mass, r-cran-robustbase, r-cran-timedate, r-cran-timeseries, r-cran-fbasics (>= 290.76), r-cran-fcopulae (>= 2100.77), xvfb, xauth, xfonts-base -Standards-Version: 3.8.4 +Build-Depends: debhelper (>= 7), r-base-dev (>= 3.0.0~), cdbs, r-cran-sn, r-cran-mass, r-cran-robustbase, r-cran-timedate, r-cran-timeseries, r-cran-fbasics (>= 290.76), r-cran-fcopulae (>= 2100.77), xvfb, xauth, xfonts-base +Standards-Version: 3.9.4 Homepage: http://www.Rmetrics.org Package: r-cran-fassets Architecture: any -Depends: ${shlibs:Depends}, r-base-core (>= 2.11.0~20100409), r-cran-sn, r-cran-mass, r-cran-robustbase, r-cran-timedate, r-cran-timeseries, r-cran-fbasics (>= 290.76), r-cran-fcopulae (>= 2100.77) +Depends: ${shlibs:Depends}, ${R:Depends}, r-cran-sn, r-cran-mass, r-cran-robustbase, r-cran-timedate, r-cran-timeseries, r-cran-fbasics (>= 290.76), r-cran-fcopulae (>= 2100.77) Suggests: r-cran-runit Description: GNU R package for financial engineering -- fAssets This package provides functions for modelling and selection of Binary files /tmp/Bvxybr48yv/fassets-2100.78/inst/DocCopying.pdf and /tmp/1NwKAyQrwV/fassets-2110.79/inst/DocCopying.pdf differ diff -Nru fassets-2100.78/man/assetsMCR.Rd fassets-2110.79/man/assetsMCR.Rd --- fassets-2100.78/man/assetsMCR.Rd 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/man/assetsMCR.Rd 2010-05-05 07:44:44.000000000 +0000 @@ -0,0 +1,151 @@ +\name{assetsMCR} + + +\alias{assetsMCR} + +\alias{covarRisk} +\alias{mcr} +\alias{mcrBeta} +\alias{riskContributions} +\alias{riskBudgets} + + +\title{Marginal Contributions to Covariance Risk} + + +\description{ + + Computes marginal contributions to covariance risk and related + measures for a portfolio of assets. + + The functions are: + + \tabular{ll}{ + \code{covarRisk} \tab Computes covariance risk, \cr + \code{mcr} \tab computes marginal contribution to covariance risk, \cr + \code{mcrBeta} \tab computes beta, the rescaled mcr to covariance risk, \cr + \code{riskConstributions} \tab computes covariance risk contributions, \cr + \code{riskBudgets} \tab computes covariance risk budgets. } + +} + + +\usage{ +covarRisk(data, weights = NULL, FUN = "cov", ...) +mcr(data, weights = NULL, FUN = "cov", ...) +mcrBeta(data, weights = NULL, FUN = "cov", ...) +riskContributions(data, weights = NULL, FUN = "cov", ...) +riskBudgets(data, weights = NULL, FUN = "cov", ...) +} + + +\arguments{ + + \item{data}{ + a multivariate 'timeSeries' object. + } + \item{weights}{ + usually a numeric vector which has the length of the number of + assets. The vector measures the weights of the individual assets. + By default \code{NULL}, then an equally weighted set of assets + is assumed. + } + \item{FUN}{ + the name of the covariance estimator function which returns the + covariance matrix. By default, the sample covariance estimator. + } + \item{\dots}{ + optional arguments to be passet to the function \code{FUN}. + } + +} + + +\value{ + + \code{covarRisk} + \cr + returns the covariance risk (standard deviation), a numeric value. + \cr + + \code{mcr} + \cr + returns the marginal contributions to covariance risk for a + portfolio of assets, a numeric value of the same length as the + number of assets. + \cr + + \code{mcrBeta} + \cr + returns the marginal contributions to beta for a + portfolio of assets, a numeric value of the same length as the + number of assets. + \cr + + \code{riskContributions} + \cr + returns the risk contributions to covariance risk for a + portfolio of assets, a numeric value of the same length as the + number of assets. + \cr + + \code{riskBudgets} + \cr + returns the risk budgets to covariance risk for a + portfolio of assets, a numeric value of the same length as the + number of assets. + +} + + +\references{ + +Goldberg L., Hayes M.Y., Menchero J., Mitra. I, (2009); + \emph{Extreme Risk Management}, + Working Paper, MSCI Barra. + +Scherer B., (2004); + \emph{Portfolio Construction and Risk Budgeting}, + Risk Books, Haymarket House. + +%Wuertz, D., Chalabi, Y., Chen W., Ellis A. (2009); +% \emph{Portfolio Optimization with R/Rmetrics}, +% Rmetrics eBook, Rmetrics Association and Finance Online, Zurich. + +} + + +\examples{ +## covarRisk - + # Covariance Risk + # Sigma = sqrt(W' COV W) + set.seed(4711) + data = assetsSim(100, 6) + covarRisk(data) + +## mcr - + # Marginal contribution to Covariance Risk + # MCR = d Sigma / d W_i + mcr(data) + +## mcrBeta - + # Marginal Beta + # beta = MCR / Sigma + mcrBeta(data) + +## riskContributions - + # Marginal Risk Contributions + # RC = Sum_i ( W_i MCR ) + riskContributions(data) + sum(riskContributions(data)) - covarRisk(data) + +## riskBudgets - + # Marginal Risk Budgets + # RB = RC / Sigma + riskBudgets(data) + sum(riskBudgets(data)) +} + + +\keyword{math} + diff -Nru fassets-2100.78/src/Makevars fassets-2110.79/src/Makevars --- fassets-2100.78/src/Makevars 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/src/Makevars 2010-05-05 07:44:43.000000000 +0000 @@ -0,0 +1 @@ +PKG_LIBS=$(BLAS_LIBS) $(FLIBS) \ No newline at end of file diff -Nru fassets-2100.78/src/ecodist.c fassets-2110.79/src/ecodist.c --- fassets-2100.78/src/ecodist.c 1970-01-01 00:00:00.000000000 +0000 +++ fassets-2110.79/src/ecodist.c 2012-09-24 07:39:48.000000000 +0000 @@ -0,0 +1,1189 @@ +#include +#include +#include /* for dgemm */ + +#define RANDIN seed_in((long *)NULL) +#define RANDOUT seed_out((long *)NULL) +#define UNIF unif_rand() +#define S_EVALUATOR + +void bootstrap(double *x, double *y, int *n, int *xlen, int *nboot, double *pboot, double *bootcor, int *rarray, int *rmat, double *xdif, double *ydif) + +{ + +int i, j, k, l; +double r; +double nsamp; +double xmean, ymean; +double xsum; +double xxsum, yysum; + +S_EVALUATOR + +/* Set random seed using Splus function */ + +RANDIN; + + +for(i = 0; i < *nboot; i++) { + +/* Set up rarray. */ + + for(j = 0; j < *n; j++) { + r = UNIF; + if(r > *pboot) + rarray[j] = 0; + else rarray[j] = 1; + } + +/* Turn rarray into a lower-triangular sampling matrix. */ +/* 1 means include, 0 means omit. */ + + l = 0; + for(j = 1; j < *n; j++) { + for(k = 0; k < j; k++) { + if(rarray[j] == 0 || rarray[k] == 0) + rmat[l] = 0; + else rmat[l] = 1; + l++; + } + } + + + nsamp = 0; + for(j = 0; j < *xlen; j++) { + nsamp += rmat[j]; + } + + +/* Calculate means for x and y. */ + + xmean = 0; + ymean = 0; + for(j = 0; j < *xlen; j++) { + if(rmat[j] == 1) { + xmean += x[j]; + ymean += y[j]; + } + } + xmean = xmean/nsamp; + ymean = ymean/nsamp; + +/* Calculate deviations for x and y. */ + + for(j = 0; j < *xlen; j++) { + if(rmat[j] == 1) { + xdif[j] = x[j] - xmean; + ydif[j] = y[j] - ymean; + } + else { + xdif[j] = 0; + ydif[j] = 0; + } + } + + xsum = 0; + xxsum = 0; + yysum = 0; + + for(j = 0; j < *xlen; j++) { + if(rmat[j] == 1) { + xsum += (xdif[j] * ydif[j]); + xxsum += (xdif[j] * xdif[j]); + yysum += (ydif[j] * ydif[j]); + } + } + + bootcor[i] = (xsum) / sqrt(xxsum * yysum); + +} + +/* Reset random seed using an Splus function. */ + +RANDOUT; + +} + + +/* DW renamed permute to permute2 */ +void permute2(double *x, double *y, int *n, int *xlen, int *nperm, + double *zstats, double *tmat, int *rarray) + +{ + +int i, k, l, m; +double cumsum; +int temp; + +S_EVALUATOR + +/* Set random seed using Splus function */ + +RANDIN; + +/* Calculate first z-statistic (unpermuted data). */ + +cumsum = 0; + +for(k = 0; k < *xlen; k++) { + cumsum += x[k] * y[k]; +} + +zstats[0] = cumsum / *xlen; + +/* Start permutation routine */ + +for(i = 1; i < *nperm; i++) { + +/* Set up rarray. */ + + for(k = 0; k < *n; k++) { + rarray[k] = k; + } + +/* Convert x to a full matrix. */ + + m = 0; + for(k = 1; k < *n; k++) { + for(l = 0; l < k; l++) { + tmat[k * *n + l] = x[m]; + tmat[l * *n + k] = x[m]; + m++; + } + } + +/* Randomize rarray using an Splus function. */ + + for(k = 0; k < (*n - 1); k++) { + l = *n - k - 1; + m = (int)((float)l * UNIF); + if(m > l) m = l; + temp = rarray[l]; + rarray[l] = rarray[m]; + rarray[m] = temp; + } + +/* Reorder x. */ + + m = 0; + for(k = 1; k < *n; k++) { + for(l = 0; l < k; l++) { + x[m] = tmat[rarray[k] * *n + rarray[l]]; + m++; + } + } + + +/* Calculate new sum of products. */ + + cumsum = 0; + + for(k = 0; k < *xlen; k++) { + cumsum += x[k] * y[k]; + + } + + zstats[i] = cumsum / *xlen; + +} + +/* Reset random seed using an Splus function. */ + +RANDOUT; + +} + + + +void permpart(double *hmat, double *bmat, double *omat, double *y, double *xcor, double *ycor, int *n, int *ncol, int *xlen, int *nperm, double *zstats, double *tmat, int *rarray) + +{ + +int i, k, l, m; +double cumsum; +double bsum; +double w1, w2; +int temp; + +S_EVALUATOR + +/* Set random seed using Splus function */ + +RANDIN; + +/* Calculate first z-statistic (unpermuted data). */ + +cumsum = 0; + +for(k = 0; k < *xlen; k++) { + cumsum += xcor[k] * ycor[k]; +} + +zstats[0] = cumsum / *xlen; + + +/* Start permutation routine */ + +for(i = 1; i < *nperm; i++) { + +/* Set up rarray. */ + + for(k = 0; k < *n; k++) { + rarray[k] = k; + } + + +/* Convert y to a full matrix. */ + + m = 0; + for(k = 1; k < *n; k++) { + for(l = 0; l < k; l++) { + tmat[k * *n + l] = y[m]; + tmat[l * *n + k] = y[m]; + m++; + } + } + +/* Randomize rarray using an Splus function. */ + + for(k = 0; k < (*n - 1); k++) { + l = *n - k - 1; + m = (int)((float)l * UNIF); + if(m > l) m = l; + temp = rarray[l]; + rarray[l] = rarray[m]; + rarray[m] = temp; + } + + +/* Reorder y. */ + + m = 0; + for(k = 1; k < *n; k++) { + for(l = 0; l < k; l++) { + y[m] = tmat[rarray[k] * *n + rarray[l]]; + m++; + } + } + + +/* Calculate residuals for y */ + +/* Calculate bmat */ + +for(k = 0; k < *ncol; k++) { + bmat[k] = 0; +} + +for(k = 0; k < *ncol; k++) { + for(l = 0; l < *xlen; l++) { + bmat[k] = bmat[k] + hmat[l * *ncol + k] * y[l]; + } +} + +/* Calculate ycor (residuals) */ + +for(k = 0; k < *xlen; k++) { + ycor[k] = 0; +} + +for(k = 0; k < *xlen; k++) { + bsum = 0; + for(l = 0; l < *ncol; l++) { + bsum = bsum + bmat[l] * omat[l * *xlen + k]; + } + ycor[k] = y[k] - bsum; +} + + +/* Standardize residuals so z = r */ + +w1 = 0; +w2 = 0; + +for(k = 0; k < *xlen; k++) { + w1 = w1 + ycor[k]; + w2 = w2 + ycor[k] * ycor[k]; +} +w1 = w1 / *xlen; +w2 = sqrt(w2 / *xlen - w1 * w1); +for(k = 0; k < *xlen; k++) { + ycor[k] = (ycor[k] - w1) / w2; +} + +/* Calculate new sum of products. */ + + cumsum = 0; + + for(k = 0; k < *xlen; k++) { + cumsum += xcor[k] * ycor[k]; + } + + zstats[i] = cumsum / *xlen; + +} + +/* Reset random seed using an Splus function. */ + +RANDOUT; + +} + + +void xbootstrap(double *x, double *y, int *n, int *xlen, int *nboot, double *pboot, double *bootcor, int *rarray, int *rmat, double *xdif, double *ydif) + +{ + +int i, j, k; +double r; +double nsamp; +double xmean, ymean; +double xsum; +double xxsum, yysum; + +S_EVALUATOR + +/* Set random seed using Splus function */ + +RANDIN; + + +for(i = 0; i < *nboot; i++) { + +/* Set up rarray. */ + + for(j = 0; j < *n; j++) { + r = UNIF; + if(r > *pboot) + rarray[j] = 0; + else rarray[j] = 1; + } + +/* Turn rarray into a square sampling matrix. */ +/* 1 means include, 0 means omit. */ + + for(j = 0; j < *xlen; j++) { + rmat[j] = 1; + } + + for(j = 0; j < *n; j++) { + for(k = 0; k <= j; k++) { + if(rarray[j] == 0 || rarray[k] == 0) { + rmat[j * *n + k] = 0; + rmat[k * *n + j] = 0; + } + } + } + + nsamp = 0; + for(j = 0; j < *xlen; j++) { + nsamp += rmat[j]; + } + + +/* Calculate means for x and y. */ + + xmean = 0; + ymean = 0; + for(j = 0; j < *xlen; j++) { + if(rmat[j] == 1) { + xmean += x[j]; + ymean += y[j]; + } + } + xmean = xmean/nsamp; + ymean = ymean/nsamp; + +/* Calculate deviations for x and y. */ + + for(j = 0; j < *xlen; j++) { + if(rmat[j] == 1) { + xdif[j] = x[j] - xmean; + ydif[j] = y[j] - ymean; + } + else { + xdif[j] = 0; + ydif[j] = 0; + } + } + + + xsum = 0; + xxsum = 0; + yysum = 0; + + for(j = 0; j < *xlen; j++) { + if(rmat[j] == 1) { + xsum += (xdif[j] * ydif[j]); + xxsum += (xdif[j] * xdif[j]); + yysum += (ydif[j] * ydif[j]); + } + } + + bootcor[i] = (xsum) / sqrt(xxsum * yysum); + +} + +/* Reset random seed using an Splus function. */ + +RANDOUT; + +} + + + +void xpermute(double *x, double *y, int *n, int *xlen, int *nperm, + double *zstats, double *tmat, int *rarray) + +{ + +int i, k, l, m; +double cumsum; +int temp; + +S_EVALUATOR + +/* Set random seed using Splus function */ + +RANDIN; + +/* Calculate first z-statistic (unpermuted data). */ + +cumsum = 0; + +for(k = 0; k < *xlen; k++) { + cumsum += x[k] * y[k]; +} + +zstats[0] = cumsum; + + +/* Start permutation routine */ + +for(i = 1; i < *nperm; i++) { + +/* Set up rarray. */ + + for(k = 0; k < *n; k++) { + rarray[k] = k; + } + +/* Randomize rarray using an Splus function. */ + + for(k = 0; k < (*n - 1); k++) { + l = *n - k - 1; + m = (int)((float)l * UNIF); + if(m > l) m = l; + temp = rarray[l]; + rarray[l] = rarray[m]; + rarray[m] = temp; + } + +/* Reorder x. */ + + for(k = 0; k < *n; k++) { + for(l = 0; l <= k; l++) { + x[k * *n + l] = tmat[rarray[k] * *n + rarray[l]]; + x[l * *n + k] = tmat[rarray[l] * *n + rarray[k]]; + } + } + + +/* Calculate new sum of products. */ + + cumsum = 0; + + for(k = 0; k < *xlen; k++) { + cumsum += x[k] * y[k]; + + } + + zstats[i] = cumsum; + +} + +/* Reset random seed using an Splus function. */ + +RANDOUT; + +} + + + + +void xpermpart(double *hmat, double *y, double *xcor, double *ycor, int *n, int *xlen, int *nperm, double *zstats, double *tmat, int *rarray) + +{ + +int i, k, l, m; +double cumsum; +int temp; + +S_EVALUATOR + +/* Set random seed using Splus function */ + +RANDIN; + +/* Calculate residuals for y */ + +for(k = 0; k < *xlen; k++) { + ycor[k] = 0; +} + + +for(k = 0; k < *xlen; k++) { + for(l = 0; l < *xlen; l++) { + ycor[k] = ycor[k] + hmat[k * *xlen + l] * y[l]; + } +} + + +/* Calculate first z-statistic (unpermuted data). */ + +cumsum = 0; + +for(k = 0; k < *xlen; k++) { + cumsum += xcor[k] * ycor[k]; +} + +zstats[0] = cumsum; + + +/* Start permutation routine */ + +for(i = 1; i < *nperm; i++) { + +/* Set up rarray. */ + + for(k = 0; k < *n; k++) { + rarray[k] = k; + } + +/* Randomize rarray using an Splus function. */ + + for(k = 0; k < (*n - 1); k++) { + l = *n - k - 1; + m = (int)((float)l * UNIF); + if(m > l) m = l; + temp = rarray[l]; + rarray[l] = rarray[m]; + rarray[m] = temp; + } + +/* Reorder y. */ + + for(k = 0; k < *n; k++) { + for(l = 0; l <= k; l++) { + y[k * *n + l] = tmat[rarray[k] * *n + rarray[l]]; + y[l * *n + k] = tmat[rarray[l] * *n + rarray[k]]; + } + } + + +/* Calculate residuals for y */ + +for(k = 0; k < *xlen; k++) { + ycor[k] = 0; +} + +for(k = 0; k < *xlen; k++) { + for(l = 0; l < *xlen; l++) { + ycor[k] = ycor[k] + hmat[k * *xlen + l] * y[l]; + } +} + + +/* Calculate new sum of products. */ + + cumsum = 0; + + for(k = 0; k < *xlen; k++) { + cumsum += xcor[k] * ycor[k]; + + } + + zstats[i] = cumsum; + +} + +/* Reset random seed using an Splus function. */ + +RANDOUT; + +} + + + +void bcdist(double *x, int *pnrow, int *pncol, double *dist) +{ +int i, j, k, l; +int nrow, ncol; +double sumi, sumj; +double minsum; + + +l = 0; +nrow = *pnrow; +ncol = *pncol; + +for(i = 0; i < (nrow - 1); i++) { + for(j = (i + 1); j < (nrow); j++) { + minsum = 0; + sumi = 0; + sumj = 0; + for(k = 0; k < ncol; k++) { + if(x[i * ncol + k] < x[j * ncol + k]) + minsum += x[i * ncol + k]; + else + minsum += x[j * ncol + k]; + sumi += x[i * ncol + k]; + sumj += x[j * ncol + k]; + } + if((sumi + sumj) == 0) + dist[l] = 0; + else + dist[l] = (1 - (2 * minsum) / (sumi + sumj)); + l++; + } +} +} + + +void weight(int *n, double *datadist, double *d1, double *d2, double *w) + +{ +int i; +double m1, m2; +double w1, w2; +double pi; + +pi = 2 * acos(0); + +for(i = 0; i < *n * *n; i++) { + if(datadist[i] != 0) { + if(d1[i] < datadist[i]) + m1 = d1[i] / datadist[i]; + else m1 = 1; + + if(d2[i] < datadist[i]) + m2 = d2[i] / datadist[i]; + else m2 = 1; + } + else { + m1 = 0; + m2 = 0; + } + + w1 = 1 - (acos(m1) + acos(m2)) / pi; + + if(datadist[i] != 0) { + m1 = d1[i] / datadist[i]; + if(m1 > 1) + m1 = 1; + + m2 = d2[i] / datadist[i]; + if(m2 > 1) + m2 = 1; + } + else { + m1 = 0; + m2 = 0; + } + + w2 = 0.75 - (acos(m1) + acos(m2)) / (2 * pi); + + if((datadist[i] * datadist[i]) >= (d1[i] * d1[i] + d2[i] * d2[i])) + w1 = 0; + if((datadist[i] * datadist[i]) < (d1[i] * d1[i] + d2[i] * d2[i])) + w2 = 0; + + w[i] = w1 + w2; +} +} + +void newpermone(double *x, int *dclass, int *n, int *xlen, int *nperm, double *zstats, double *tmat, int *rarray) + +{ + +int i, k, l, m; +double cumsum; +int temp; + +S_EVALUATOR + +/* Set random seed using Splus function */ + +RANDIN; + +/* Calculate first z-statistic (unpermuted data). */ + +cumsum = 0; + +for(k = 0; k < *xlen; k++) { + if(dclass[k] == 0) { + cumsum += x[k]; + } +} + +zstats[0] = cumsum; + +/* Start permutation routine */ + +for(i = 1; i < *nperm; i++) { + +/* Set up rarray. */ + + for(k = 0; k < *n; k++) { + rarray[k] = k; + } + +/* Convert x to a full matrix. */ + + m = 0; + for(k = 1; k < *n; k++) { + for(l = 0; l < k; l++) { + tmat[k * *n + l] = x[m]; + tmat[l * *n + k] = x[m]; + m++; + } + } + +/* Randomize rarray using an Splus function. */ + + for(k = 0; k < (*n - 1); k++) { + l = *n - k - 1; + m = (int)((float)l * UNIF); + if(m > l) m = l; + temp = rarray[l]; + rarray[l] = rarray[m]; + rarray[m] = temp; + } + +/* Reorder x. */ + + m = 0; + for(k = 1; k < *n; k++) { + for(l = 0; l < k; l++) { + x[m] = tmat[rarray[k] * *n + rarray[l]]; + m++; + } + } + + +/* Calculate new sum of products. */ + + cumsum = 0; + + for(k = 0; k < *xlen; k++) { + if(dclass[k] == 0) { + cumsum += x[k]; + } + } + + zstats[i] = cumsum; + +} + +/* Reset random seed using an Splus function. */ + +RANDOUT; + +} + + + + +void newpermtwo(double *x, double *y, int *n, int *xlen, int *nperm, double *zstats, double *tmat, int *rarray) + +{ + +int i, k, l, m; +double cumsum; +int temp; +float naval = -9999; + +S_EVALUATOR + +/* Set random seed using Splus function */ + +RANDIN; + +/* Calculate first z-statistic (unpermuted data). */ + +cumsum = 0; + +for(k = 0; k < *xlen; k++) { + if(x[k] != naval) { + cumsum += x[k] * y[k]; + } +} + +zstats[0] = cumsum; + +/* Start permutation routine */ + +for(i = 1; i < *nperm; i++) { + +/* Set up rarray. */ + + for(k = 0; k < *n; k++) { + rarray[k] = k; + } + +/* Convert x to a full matrix. */ + + m = 0; + for(k = 1; k < *n; k++) { + for(l = 0; l < k; l++) { + tmat[k * *n + l] = x[m]; + tmat[l * *n + k] = x[m]; + m++; + } + } + +/* Randomize rarray using an Splus function. */ + + for(k = 0; k < (*n - 1); k++) { + l = *n - k - 1; + m = (int)((float)l * UNIF); + if(m > l) m = l; + temp = rarray[l]; + rarray[l] = rarray[m]; + rarray[m] = temp; + } + +/* Reorder x. */ + + m = 0; + for(k = 1; k < *n; k++) { + for(l = 0; l < k; l++) { + x[m] = tmat[rarray[k] * *n + rarray[l]]; + m++; + } + } + + +/* Calculate new sum of products. */ + + cumsum = 0; + + for(k = 0; k < *xlen; k++) { + if(x[k] != naval) { + cumsum += x[k] * y[k]; + } + } + + zstats[i] = cumsum; + +} + +/* Reset random seed using an Splus function. */ + +RANDOUT; + +} + + + + +void psum(double *x, int *pnrow, int *pncol, double *dist) +{ +int row1, row2, col1; +int nrow, ncol; +int l; +double thisval, thatval; + +l = 0; +nrow = *pnrow; +ncol = *pncol; + + for(col1 = 0; col1 < ncol; col1++) { + for(row1 = 0; row1 < nrow; row1++) { + thatval = x[row1 * ncol + col1]; + for(row2 = 0; row2 < nrow; row2++) { + thisval = x[row2 * ncol + col1]; + dist[l] = thisval + thatval; + l++; + } + } +} + +} + + +void pdiff(double *x, int *pnrow, int *pncol, double *dist) +{ +int row1, row2, col1; +int nrow, ncol; +int l; +double thisval, thatval; + +l = 0; +nrow = *pnrow; +ncol = *pncol; + + for(col1 = 0; col1 < ncol; col1++) { + for(row1 = 0; row1 < nrow; row1++) { + thatval = x[row1 * ncol + col1]; + for(row2 = 0; row2 < nrow; row2++) { + thisval = x[row2 * ncol + col1]; + dist[l] = thisval - thatval; + l++; + } + } +} + +} + + +void jpres(double *x, int *pnrow, int *pncol, double *dist) +{ +int row1, row2, col1; +int nrow, ncol; +int l; +double thisval, thatval; + +l = 0; +nrow = *pnrow; +ncol = *pncol; + + for(col1 = 0; col1 < ncol; col1++) { + for(row1 = 0; row1 < nrow; row1++) { + thatval = x[row1 * ncol + col1]; + for(row2 = 0; row2 < nrow; row2++) { + thisval = x[row2 * ncol + col1]; + if((thisval > 0) & (thatval > 0)) { + dist[l] = 1; + } + else { + dist[l] = 0; + } + l++; + } + } +} + +} + + +void jabs(double *x, int *pnrow, int *pncol, double *dist) +{ +int row1, row2, col1; +int nrow, ncol; +int l; +double thisval, thatval; + +l = 0; +nrow = *pnrow; +ncol = *pncol; + + for(col1 = 0; col1 < ncol; col1++) { + for(row1 = 0; row1 < nrow; row1++) { + thatval = x[row1 * ncol + col1]; + for(row2 = 0; row2 < nrow; row2++) { + thisval = x[row2 * ncol + col1]; + if((thisval == 0) & (thatval == 0)) { + dist[l] = 1; + } + else { + dist[l] = 0; + } + l++; + } + } +} + +} + + +void jfirst(double *x, int *pnrow, int *pncol, double *dist) +{ +int row1, row2, col1; +int nrow, ncol; +int l; +double thisval, thatval; + +l = 0; +nrow = *pnrow; +ncol = *pncol; + + for(col1 = 0; col1 < ncol; col1++) { + for(row1 = 0; row1 < nrow; row1++) { + thatval = x[row1 * ncol + col1]; + for(row2 = 0; row2 < nrow; row2++) { + thisval = x[row2 * ncol + col1]; + if((thisval > 0) & (thatval == 0)) { + dist[l] = 1; + } + else { + dist[l] = 0; + } + l++; + } + } +} + +} + + +void jsec(double *x, int *pnrow, int *pncol, double *dist) +{ +int row1, row2, col1; +int nrow, ncol; +int l; +double thisval, thatval; + +l = 0; +nrow = *pnrow; +ncol = *pncol; + + for(col1 = 0; col1 < ncol; col1++) { + for(row1 = 0; row1 < nrow; row1++) { + thatval = x[row1 * ncol + col1]; + for(row2 = 0; row2 < nrow; row2++) { + thisval = x[row2 * ncol + col1]; + if((thisval == 0) & (thatval > 0)) { + dist[l] = 1; + } + else { + dist[l] = 0; + } + l++; + } + } +} + +} + + +void mrmperm(double *x, double *y, int *p, int *nd, int *n, int *nperm, double *r2all, double *ball, double *fall, double *tmat, int *rarray, double *XX, double *XY, double *YY, double *b) + +{ + +int i, k, l; +int m; +int temp; +double SSE=0.0, SSTO=0.0, SSR=0.0; +double r2=0, f=0; +double btemp=0.0; +int bcount = 0; +char *transt = "T", *transn = "N"; +double one = 1.0, zero = 0.0; +int onei = 1; + +S_EVALUATOR + +/* Set random seed using Splus function */ + +RANDIN; + +/* Start permutation routine */ +for(i = 0; i < *nperm; i++) { + +/* first do the unpermuted values */ + +/* F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one, + x, &nrx, y, &nry, &zero, z, &ncx); */ + +/* take crossproduct t(X) %*% Y - WORKS */ + F77_CALL(dgemm)(transt, transn, + p, &onei, nd, + &one, x, nd, y, nd, + &zero, XY, p); + +/* take crossproduct t(Y) %*% (Y) - WORKS */ + F77_CALL(dgemm)(transt, transn, + &onei, &onei, nd, + &one, y, nd, y, nd, + &zero, YY, &onei); + +/* calculate regression coefficients XX %*% XY - WORKS */ + F77_CALL(dgemm)(transn, transn, + p, &onei, p, + &one, XX, p, XY, p, + &zero, b, p); + +/* calculate regression components - WORKS */ + F77_CALL(dgemm)(transt, transn, + &onei, &onei, p, + &one, b, p, XY, p, + &zero, &btemp, &onei); + +/* SSE - WORKS */ + SSE = YY[0] - btemp; + +/* SSTO - WORKS */ + SSTO = 0; + for(k = 0; k < *nd; k++) { + SSTO = SSTO + y[k]; + } + SSTO = YY[0] - (SSTO * SSTO) / *nd; + + SSR = SSTO - SSE; + + /* calculate R2 - WORKS */ + r2 = 1 - SSE / SSTO; + + /* calculate F* - WORKS */ + f = (SSR / (*p - 1)) / (SSE / (*nd - *p)); + + r2all[i] = r2; + fall[i] = f; + + /* calculate pseudo-t for regression coefficients - WORKS*/ + /* b / sqrt(1 - R2) */ + for(k=0; k<*p; k++) { + ball[bcount] = b[k] / sqrt(1 - r2); + bcount++; + } + + +/* permute Y */ +/* Set up rarray. */ + + for(k = 0; k < *n; k++) { + rarray[k] = k; + } + +/* Convert y to a full matrix. */ + + m = 0; + for(k = 1; k < *n; k++) { + for(l = 0; l < k; l++) { + tmat[k * *n + l] = y[m]; + tmat[l * *n + k] = y[m]; + m++; + } + } + +/* Randomize rarray using an Splus function. */ + + for(k = 0; k < (*n - 1); k++) { + l = *n - k - 1; + m = (int)((float)l * UNIF); + if(m > l) m = l; + temp = rarray[l]; + rarray[l] = rarray[m]; + rarray[m] = temp; + } + +/* Reorder y. */ + + m = 0; + for(k = 1; k < *n; k++) { + for(l = 0; l < k; l++) { + y[m] = tmat[rarray[k] * *n + rarray[l]]; + m++; + } + } + +} + +/* Reset random seed using an Splus function. */ + +RANDOUT; + +} + +