[Returnanalytics-commits] r3682 - in pkg/Meucci: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jun 20 10:11:15 CEST 2015
Author: xavierv
Date: 2015-06-20 10:11:15 +0200 (Sat, 20 Jun 2015)
New Revision: 3682
Modified:
pkg/Meucci/NAMESPACE
pkg/Meucci/R/DetectOutliersviaMVE.R
pkg/Meucci/R/EntropyProg.R
pkg/Meucci/R/FullyFlexibleBayesNets.R
pkg/Meucci/R/HermiteGrid.R
pkg/Meucci/R/data.R
pkg/Meucci/demo/DetectOutliersviaMVE.R
pkg/Meucci/demo/FullFlexProbs.R
pkg/Meucci/demo/FullyFlexibleBayesNets.R
pkg/Meucci/demo/FullyIntegratedLiquidityAndMarketRisk.R
pkg/Meucci/demo/HermiteGrid_CVaR_Recursion.R
pkg/Meucci/demo/HermiteGrid_CaseStudy.R
pkg/Meucci/demo/HermiteGrid_demo.R
pkg/Meucci/demo/RobustBayesianAllocation.R
pkg/Meucci/man/ComputeMVE.Rd
pkg/Meucci/man/ComputeMoments.Rd
pkg/Meucci/man/CondProbViews.Rd
pkg/Meucci/man/DetectOutliersViaMVE.Rd
pkg/Meucci/man/Equities.Rd
pkg/Meucci/man/NoisyObservations.Rd
pkg/Meucci/man/RejectOutlier.Rd
pkg/Meucci/man/StockSeries.Rd
pkg/Meucci/man/TimeSeries.Rd
pkg/Meucci/man/Tweak.Rd
pkg/Meucci/man/UsSwapRates.Rd
pkg/Meucci/man/bondAttribution.Rd
pkg/Meucci/man/butterfliesAnalytics.Rd
pkg/Meucci/man/covNRets.Rd
pkg/Meucci/man/db.Rd
pkg/Meucci/man/dbFFP.Rd
pkg/Meucci/man/db_FX.Rd
pkg/Meucci/man/derivatives.Rd
pkg/Meucci/man/fILMR.Rd
pkg/Meucci/man/factorsDistribution.Rd
pkg/Meucci/man/fixedIncome.Rd
pkg/Meucci/man/freaqEst.Rd
pkg/Meucci/man/gaussHermiteMesh.Rd
pkg/Meucci/man/hermitePolynomial.Rd
pkg/Meucci/man/highYieldIndices.Rd
pkg/Meucci/man/implVol.Rd
pkg/Meucci/man/integrateSubIntervals.Rd
pkg/Meucci/man/kernelbw.Rd
pkg/Meucci/man/kernelcdf.Rd
pkg/Meucci/man/kernelinv.Rd
pkg/Meucci/man/kernelpdf.Rd
pkg/Meucci/man/linearModel.Rd
pkg/Meucci/man/normalizeProb.Rd
pkg/Meucci/man/private_fun.Rd
pkg/Meucci/man/returnsDistribution.Rd
pkg/Meucci/man/sectorsSnP500.Rd
pkg/Meucci/man/sectorsTS.Rd
pkg/Meucci/man/securitiesIndustryClassification.Rd
pkg/Meucci/man/securitiesTS.Rd
pkg/Meucci/man/subIntervals.Rd
pkg/Meucci/man/swap2y4y.Rd
pkg/Meucci/man/swapParRates.Rd
pkg/Meucci/man/swaps.Rd
Log:
- Formatted demo scripts and its functions up to HermiteGrid_demo
Modified: pkg/Meucci/NAMESPACE
===================================================================
--- pkg/Meucci/NAMESPACE 2015-06-18 22:35:26 UTC (rev 3681)
+++ pkg/Meucci/NAMESPACE 2015-06-20 08:11:15 UTC (rev 3682)
@@ -43,6 +43,7 @@
export(MvnRnd)
export(NoisyObservations)
export(NormalCopulaPdf)
+export(PHist)
export(PanicCopula)
export(PartialConfidencePosterior)
export(PerformIidAnalysis)
Modified: pkg/Meucci/R/DetectOutliersviaMVE.R
===================================================================
--- pkg/Meucci/R/DetectOutliersviaMVE.R 2015-06-18 22:35:26 UTC (rev 3681)
+++ pkg/Meucci/R/DetectOutliersviaMVE.R 2015-06-20 08:11:15 UTC (rev 3682)
@@ -1,16 +1,18 @@
-#' Finds the "worst" outlier in a multivariate time series
+#' @title Finds the "worst" outlier in a multivariate time series
#'
-#' Finds the "worst" outlier in a multivariate time-series
+#' @description Finds the "worst" outlier in a multivariate time series.
+#'
#' We aim at finding the the observation x_t such that if we remove it from the
-#' set {x_1 ... x_t } the determinant of the resulting sample covariance is reduced the most
-#' This means that by dropping that observation the location-dispersion ellipsoid defined
-#' by the sample mean and covariance shrinks the most
-#' See Sec. 4.6.1 of "Risk and Asset Allocation" - Springer (2005), by A. Meucci
-#' for the theory and the routine implemented below
+#' set {x_1 ... x_t } the determinant of the resulting sample covariance is
+#' reduced the most. This means that by dropping that observation the
+#' location-dispersion ellipsoid defined by the sample mean and covariance
+#' shrinks the most. See Sec. 4.6.1 of "Risk and Asset Allocation" - Springer
+#' (2005), by A. Meucci for the theory and the routine implemented below
#'
#' @param sample a vector containing the input dataset with outliers
#'
-#' @return rejected a numeric indicating which observation in the index to reject
+#' @return rejected a numeric indicating which observation in the index to
+#' reject
#'
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
#'
@@ -18,29 +20,35 @@
#' \url{http://www.symmys.com}
#' See Meucci script for "RejectOutlier.m"
#' @export
-RejectOutlier = function( sample )
-{
- library( matlab )
-
+
+RejectOutlier <- function(sample) {
+ library(matlab)
# parameter checks
- if ( ncol( sample ) > nrow( sample ) ) { stop("The number of assets must be greater than number of observations (otherwise system is singular)") }
-
+ if (ncol(sample) > nrow(sample)) {
+ stop("The number of assets must be greater than number of observations
+ (otherwise system is singular)")
+ }
+
# initialize parameters
- T = nrow( sample )
- m = matrix( colMeans( sample ) , nrow = ncol( sample ) )
- U = sample - ones( T , 1 ) %*% t( m )
-
+ T <- nrow(sample)
+ m <- matrix(colMeans(sample), nrow = ncol(sample))
+ U <- sample - ones(T, 1) %*% t(m)
+
# measure lambdas associated with each observation
- lambdas = diag( U %*% solve( t(U) %*% U ) %*% t(U) ) # result is from Poston, Wegman, Priebe and Solka (1997)
- # the various lambdas denote the t-th element of the diagonal of the information matrix
+
+ # result is from Poston, Wegman, Priebe and Solka (1997)
+ lambdas <- diag(U %*% solve(t(U) %*% U) %*% t(U))
+ # the various lambdas denote the t-th element of the diagonal of the
+ # information matrix
# the farthest outlier corresponds to the highest value of lambda
- a = max( lambdas )
- rejected = which.max( lambdas ) # identify the farthest outlier
-
- if ( max(lambdas > 1) == 1 ) { stop( "Lambdas cannot be greater than 1")}
- if ( max(lambdas < 0) == 1 ) { stop( "Lambdas cannot be less than 0")}
-
- return( list( rejected = rejected , lambdas = lambdas ) )
+ rejected <- which.max(lambdas) # identify the farthest outlier
+
+ if (max(lambdas > 1) == 1)
+ stop("Lambdas cannot be greater than 1")
+ if (max(lambdas < 0) == 1)
+ stop("Lambdas cannot be less than 0")
+
+ return(list(rejected = rejected, lambdas = lambdas))
}
#' Compute the minimum volume ellipsoid for a given (multi-variate) time-series
@@ -49,7 +57,7 @@
#'
#' via the expectations-minimization algorithm
#'
-#' \deqn{ w_{t} = \frac{1}{T} , t = 1,...,T
+#' \deqn{ w_{t} = \frac{1}{T}, t = 1,...,T
#' \\ m \equiv \frac{1}{ \sum_{s=1}^T w_{s} } \sum_{t=1}^T w_{t} x_{t}
#' \\ S \equiv \sum_{t=1}^T w_{t} \big(x_{t} - m\big) \big(x_{t} - m\big)'
#' \\ Ma_{t}^{2} \equiv \big(x-m\big)' S^{-1} \big(x-m\big), t=1,...,T
@@ -60,171 +68,199 @@
#' The location and scatter parameters that define the ellipsoid are
#' multivariate high-breakdown estimators of location and scatter
#'
-#' @param data a matrix time-series of data. Each row is a observation (date). Each column is an asset
-#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
+#' @param data a matrix time-series of data. Each row is a observation
+#' (date). Each column is an asset
+#'
#' @return list a list with
-#' MVE_Location a numeric with the location parameter of minimum volume ellipsoid
-#' MVE_Dispersion a numeric with the covariance matrix of the minimum volume ellipsoid
+#' MVE_Location a numeric with the location parameter of minimum volume
+#' ellipsoid
+#' MVE_Dispersion a numeric with the covariance matrix of the minimum
+#' volume ellipsoid
#'
#' @references
-#' \url{http://www.symmys.com/sites/default/files/Risk\%20and\%20Asset\%20Allocation\%20-\%20Springer\%20Quantitative\%20Finance\%20-\%20Estimation.pdf}
#' See Sec. 4.6.1 of "Risk and Asset Allocation" - Springer (2005), by A. Meucci
#' for the theory and the routine implemented below
#' See Meucci script for "ComputeMVE.m"
+#'
+#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
+#'
#' @export
-#'
-ComputeMVE = function ( data )
-{
- library( matlab )
- NumObservations = nrow( data )
- Ones = ones( NumObservations , 1 )
- det_S_New = 0
-
+
+ComputeMVE <- function (data) {
+ library(matlab)
+ NumObservations <- nrow(data)
+ Ones <- ones(NumObservations, 1)
+ det_S_New <- 0
+
# step 0: initialize the relative weights
- w = 1 / NumObservations * Ones
-
+ w <- 1 / NumObservations * Ones
+
# step 1: compute the location parameter (m), and the scatter matrix (S)
- m = matrix( colMeans( data ) , nrow = ncol( data ) )
- S = cov( data ) # notice the weights in the scatter matrix are not normalized
-
- keeploop = TRUE
- while ( keeploop == TRUE )
- {
- Mahalanobis = matrix( , ncol = 0 , nrow = 1 )
- for ( t in 1:NumObservations )
- {
+ m <- matrix(colMeans(data), nrow = ncol(data))
+ S <- cov(data) # notice the weights in the scatter matrix are not normalized
+
+ keeploop <- TRUE
+ while (keeploop == TRUE) {
+ Mahalanobis <- matrix(, ncol = 0, nrow = 1)
+ for (t in 1:NumObservations) {
# cycle thru each observation...
- x_t = t( data[ t , , drop = FALSE ] )
-
+ x_t <- t(data[t, , drop = FALSE])
+
# ...and calculate the square Mahalanobis distances
- Mahalanobis = cbind( Mahalanobis , t((x_t - m)) %*% solve(S) %*% ( x_t - m ) )
- }
- # Step 3: update the weights if the Mahalanobis distance squared > 1 as follows...
- update = matlab:::find(Mahalanobis > 1 )
- w[ update ] = w[ update ] * t( Mahalanobis[ update ] )
+ Mahalanobis <- cbind(Mahalanobis, t((x_t - m)) %*% solve(S) %*% (x_t - m))
+ }
+ # Step 3: update the weights if Mahalanobis distance squared > 1 as follows
+ update <- matlab:::find(Mahalanobis > 1)
+ w[update] <- w[update] * t(Mahalanobis[update])
# ... otherwise leave the weight unchanged
-
- # Step 4: If convergence is reached, stop and define mu and Sigma, otherwise go to Step 1
- m = t( data ) %*% w / sum (w )
- S = t(( data - Ones %*% t( m ) )) %*% diag( as.vector( w ) ) %*% ( data - Ones %*% t(m) )
-
- det_S_Old = det_S_New
- det_S_New = base:::det(S)
-
- exitFlag = ( ( det_S_Old / det_S_New ) < .99999 ) # loop exits when evaluates to FALSE
-
- if ( det_S_New == 0 ) { exitFlag = TRUE }
-
- if ( !is.logical( exitFlag ) ) { browser() }
-
- keeploop = ( exitFlag )
- }
-
- MVE_Location = m
- MVE_Dispersion = S
-
- return( list( MVE_Location = MVE_Location , MVE_Dispersion = MVE_Dispersion ) )
+
+ # Step 4: If convergence is reached, stop and define mu and Sigma,
+ # otherwise go to Step 1
+ m <- t(data) %*% w / sum (w)
+ S <- t((data - Ones %*% t(m))) %*% diag(as.vector(w)) %*% (data -
+ Ones %*% t(m))
+
+ det_S_Old <- det_S_New
+ det_S_New <- base:::det(S)
+
+ # loop exits when evaluates to FALSE
+ exitFlag <- ((det_S_Old / det_S_New) < .99999)
+
+ if (det_S_New == 0)
+ exitFlag <- TRUE
+
+ if (!is.logical(exitFlag))
+ browser()
+
+ keeploop <- (exitFlag)
+ }
+
+ MVE_Location <- m
+ MVE_Dispersion <- S
+
+ return(list(MVE_Location = MVE_Location, MVE_Dispersion = MVE_Dispersion))
}
-#' Use the minimum volume ellipsoid to detect outliers
+#' @title Use the minimum volume ellipsoid to detect outliers
#'
-#' See Sec. 4.6.1 of "Risk and Asset Allocation" - Springer (2005), by A. Meucci
-#' for the theory and the routine implemented below
+#' @description Use the minimum volume ellipsoid to detect outliers. See Sec.
+#' 4.6.1 of "Risk and Asset Allocation" - Springer (2005), by A. Meucci for the
+#' theory and the routine implemented below.
#'
-#' @param corruptSample a matrix of returns with outlier data. Rows are observations, columns are assets.
+#' @param corruptSample a matrix of returns with outlier data. Rows are
+#' observations, columns are assets.
#'
#' @return a list containing:
-#' plotdata a matrix of data used to plot minimum volume ellipsoid as a function of its length
-#' cutofflist an ordering of observations with the highest Mahalanobis distance (i.e. ordering of outliers by their index )#'
-#' numOutliers returns the number of outliers based on the slope of the minimum volume ellipsoid as a function of sample data
+#' plotdata a matrix of data used to plot minimum volume ellipsoid as a
+#' function of its length
+#' cutofflist an ordering of observations with the highest Mahalanobis
+#' distance (i.e. ordering of outliers by their index)
+#' numOutliers returns the number of outliers based on the slope of the
+#' minimum volume ellipsoid as a function of sample data
#'
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
#' @references
#' \url{http://www.symmys.com}
#' See Meucci script for "S_HighBreakdownMVE.m"
#' @export
-DetectOutliersViaMVE = function( corruptSample )
-{
- library( matlab )
-
+
+DetectOutliersViaMVE <- function(corruptSample) {
+
# parameter checks
- if ( ncol( corruptSample ) > nrow( corruptSample ) ) { stop("The number of assets must be greater than number of observations (otherwise system is singular)") }
-
+ if (ncol(corruptSample) > nrow(corruptSample)) {
+ stop("The number of assets must be greater than number of observations
+ (otherwise system is singular)")
+ }
+
# initialize variables
- T = nrow( corruptSample )
- index = seq( from = 1 , to = T , 1 )
- vol_MVE = sample_Length = matrix( , nrow = 1 , ncol = 0 )
- store = rep( list() , length = ceil( T / 2 ) )
-
- lambdas = RejectOutlier( corruptSample )$lambdas
- cutofflist = order( lambdas , decreasing = TRUE )
-
+ T <- nrow(corruptSample)
+ index <- seq(from = 1, to = T, 1)
+ vol_MVE <- sample_Length <- matrix(, nrow = 1, ncol = 0)
+ store <- rep(list(), length = ceil(T / 2))
+
+ lambdas <- RejectOutlier(corruptSample)$lambdas
+ cutofflist <- order(lambdas, decreasing = TRUE)
+
# compute high-breakdown estimates
- for ( j in 1:ceil( T / 2 ) ) # stop when the number of data left is less than half the original sample, otherwise repeat
- {
- # step 1: compute MVE location and dispersion
- MVE = ComputeMVE( corruptSample )
-
+
+ # stop when the number of data left is less than half the original sample,
+ # otherwise repeat
+ for (j in 1:ceil(T / 2)) {
+ # step 1: compute MVE location and dispersion
+ MVE <- ComputeMVE(corruptSample)
+
# store results
- store[[j]]$MVE_Location = MVE$MVE_Location
- store[[j]]$MVE_Dispersion = MVE$MVE_Dispersion
-
- sample_Length = cbind( sample_Length , nrow( corruptSample ) )
- vol_MVE = cbind( vol_MVE , sqrt( det( MVE$MVE_Dispersion ) ) )
-
- store[[j]]$index = index
-
- # step 2: find the farthest outlier among the data and remove the observation (one at a time)
- rejected = RejectOutlier( corruptSample )$rejected
- corruptSample = corruptSample[ -rejected , ] # drop outlier rows from sample
- index = index[ -rejected ]
+ store[[j]]$MVE_Location <- MVE$MVE_Location
+ store[[j]]$MVE_Dispersion <- MVE$MVE_Dispersion
+
+ sample_Length <- cbind(sample_Length, nrow(corruptSample))
+ vol_MVE <- cbind(vol_MVE, sqrt(det(MVE$MVE_Dispersion)))
+
+ store[[j]]$index <- index
+
+ # step 2: find the farthest outlier among the data and remove the
+ # observation (one at a time)
+ rejected <- RejectOutlier(corruptSample)$rejected
+ corruptSample <- corruptSample[-rejected, ] # drop outlier rows from sample
+ index <- index[-rejected]
}
-
- # plot volume of ellipsoid as function of sample length
+
+ # plot volume of ellipsoid as function of sample length
# shows an abrupt jump when the first outlier is added to the data
- # The respective sample covariance is the minimum covariance determinant and the respective ellipsoid
- plot( sample_Length , vol_MVE , type = "l" , main = "Outlier detection" , ylab = "volume of Min volume ellipsoid" , xlab = "sample length" )
-
+ # The respective sample covariance is the minimum covariance determinant and
+ # the respective ellipsoid
+ plot(sample_Length, vol_MVE, type = "l", main = "Outlier detection",
+ ylab = "volume of Min volume ellipsoid", xlab = "sample length")
+
# return the minimum volume ellipsoid at each sample length
- result = rbind( sample_Length , vol_MVE )
- rownames( result ) = c("index" , "volumeMVE" )
-
- # TODO: measure the slope of result$volumeMVE. When slope doubles that marks the first outlier
- numOutliers = NULL
-
- return( list( plotdata = t( result ) , cutofflist = cutofflist , numOutliers = numOutliers ) )
+ result <- rbind(sample_Length, vol_MVE)
+ rownames(result) = c("index", "volumeMVE")
+
+ # TODO: measure the slope of result$volumeMVE.
+ # When slope doubles that marks the first outlier
+ numOutliers <- NULL
+
+ return(list(plotdata = t(result), cutofflist = cutofflist,
+ numOutliers = numOutliers))
}
-#' Generate observations from a two asset covariance matrix and add outliers
+#' @title Generate observations from a covariance matrix and add outliers
#'
-#' Generate observations from a covariance matrix and add outliers
+#' @description Generate observations from a two asset covariance matrix and add outliers
#'
-#' @param numGoodSamples number of observations drawn from the covariance matrix
+#' @param numGoodSamples number of observations drawn from the
+#' covariance matrix
#' @param numOutliers number of outliers added to sample
-#' @param covarianceMatrix the covariance matrix for the asset returns from which good samples will be drawn
-#' @param shuffle a boolean suggesting whether order of the twos should be shuffled
+#' @param covarianceMatrix the covariance matrix for the asset returns
+#' from which good samples will be drawn
+#' @param shuffle a boolean suggesting whether order of the twos
+#' should be shuffled
#'
-#' @return sample a matrix of returns consisting of good and bad sample. Rows are observations, columns are the assets.
+#' @return sample a matrix of returns consisting of good and bad sample.
+#' Rows are observations, columns are the assets.
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
#' @export
-NoisyObservations = function( numGoodSamples , numOutliers , covarianceMatrix , shuffle = FALSE )
-{
- mu = matrix( rep( 0 , nrow( covarianceMatrix ) ) )
- T = numGoodSamples
-
+
+NoisyObservations <- function(numGoodSamples, numOutliers, covarianceMatrix,
+ shuffle = FALSE) {
+ mu <- matrix(rep(0, nrow(covarianceMatrix)))
+ T <- numGoodSamples
+ require("MASS")
# "good" observations
- library( MASS )
- sample = mvrnorm( n = T , mu = mu , Sigma = covarianceMatrix ) # draw T samples from mvrnorm distribution
-
+
+ # draw T samples from mvrnorm distribution
+ sample <- mvrnorm(n = T, mu = mu, Sigma = covarianceMatrix)
# generate outliers
- outliers = 10 * matrix( runif( numOutliers * nrow( covarianceMatrix ) ) , nrow = numOutliers , ncol = ncol( covarianceMatrix ) )
-
+ outliers <- 10 * matrix(runif(numOutliers * nrow(covarianceMatrix)),
+ nrow = numOutliers, ncol = ncol(covarianceMatrix))
+
# add "bad" observation(s)
- corruptSample = rbind ( sample , outliers )
-
+ corruptSample <- rbind (sample, outliers)
+
# shuffle rows
- if ( shuffle == TRUE ) { corruptSample = corruptSample[ sample( nrow( corruptSample ) ) , ] }
-
- return( corruptSample = corruptSample )
+ if (shuffle == TRUE)
+ corruptSample <- corruptSample[sample(nrow(corruptSample)), ]
+
+ return(corruptSample = corruptSample)
}
Modified: pkg/Meucci/R/EntropyProg.R
===================================================================
--- pkg/Meucci/R/EntropyProg.R 2015-06-18 22:35:26 UTC (rev 3681)
+++ pkg/Meucci/R/EntropyProg.R 2015-06-20 08:11:15 UTC (rev 3682)
@@ -311,6 +311,8 @@
#' See Meucci script pHist.m used for plotting
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com} and Xavier Valls
#' \email{xaviervallspla@@gmail.com}
+#'
+#' @export
PHist <- function(X, p, nBins, freq = FALSE) {
Modified: pkg/Meucci/R/FullyFlexibleBayesNets.R
===================================================================
--- pkg/Meucci/R/FullyFlexibleBayesNets.R 2015-06-18 22:35:26 UTC (rev 3681)
+++ pkg/Meucci/R/FullyFlexibleBayesNets.R 2015-06-20 08:11:15 UTC (rev 3682)
@@ -1,9 +1,10 @@
#' Input conditional views
#'
#' statement: View(k).Who (e.g. [1 3])= View(k).Equal (e.g. {[2 3] [1 3 5]})
-#' optional conditional statement: View(k).Cond_Who (e.g. [2])= View(k).Cond_Equal (e.g. {[1]})
-#' amount of stress is quantified as Prob(statement) <= View(k).v if View(k).sgn = 1;
-#' Prob(statement) >= View(k).v if View(k).sgn = -1;
+#' optional conditional statement: View(k).Cond_Who (e.g. [2]) =
+#' View(k).Cond_Equal (e.g. {[1]})
+#' amount of stress quantified as Prob(statement)<=View(k).v if View(k).sgn = 1
+#' Prob(statement)>=View(k).v if View(k).sgn = -1
#'
#' confidence in stress is quantified in View(k).c in (0,1)
#'
@@ -15,145 +16,169 @@
#' @return g TBD
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
#' @export
-CondProbViews = function( View , X )
-{
- # initialize parameters
- A = matrix( , nrow = 0 , ncol = nrow( X ) )
- b = g = matrix( , nrow = 0 , ncol = 1 )
-
- # for each view...
- for ( k in 1:length( View ) ) {
- I_mrg = ( X[ , 1] < Inf )
-
- for ( s in 1:length( View[[k]]$Who ) )
- {
- Who = View[[k]]$Who[s]
- Or_Targets = View[[k]]$Equal[[s]]
- I_mrg_or = ( X[ , Who] > Inf )
- for ( i in 1:length( Or_Targets ) ) { I_mrg_or = I_mrg_or | ( X[ , Who ] == Or_Targets[i] ) } # element-wise logical OR
- I_mrg = I_mrg & I_mrg_or # element-wise logical AND
+
+CondProbViews <- function(View, X) {
+ # initialize parameters
+ A <- matrix(, nrow = 0, ncol = nrow(X))
+ b <- g <- matrix(, nrow = 0, ncol = 1)
+
+ # for each view...
+ for (k in 1:length(View)) {
+ I_mrg <- (X[, 1] < Inf)
+
+ for (s in 1:length(View[[k]]$Who)) {
+ Who <- View[[k]]$Who[s]
+ Or_Targets <- View[[k]]$Equal[[s]]
+ I_mrg_or <- (X[, Who] > Inf)
+ for (i in 1:length(Or_Targets)) {
+ I_mrg_or <- I_mrg_or | (X[, Who] == Or_Targets[i])
+ } # element-wise logical OR
+ I_mrg <- I_mrg & I_mrg_or # element-wise logical AND
}
-
- I_cnd = ( X[ , 1 ] < Inf )
-
- if ( length( View[[k]]$Cond_Who ) != 0) # If length of Cond_Who is zero, skip
- {
- for ( s in 1:length( View[[k]]$Cond_Who ) )
- {
- Who = View[[k]]$Cond_Who[s]
- Or_Targets = View[[k]]$Cond_Equal[[s]]
- I_cnd_or = ( X[ , Who ] > Inf )
- for ( i in 1:length( Or_Targets ) ) { I_cnd_or = I_cnd_or | X[ ,Who ] == Or_Targets[i] }
- I_cnd = I_cnd & I_cnd_or
+
+ I_cnd <- (X[, 1] < Inf)
+ # If length of Cond_Who is zero, skip
+ if (length(View[[k]]$Cond_Who) != 0) {
+ for (s in 1:length(View[[k]]$Cond_Who)) {
+ Who <- View[[k]]$Cond_Who[s]
+ Or_Targets <- View[[k]]$Cond_Equal[[s]]
+ I_cnd_or <- (X[, Who] > Inf)
+ for (i in 1:length(Or_Targets)) {
+ I_cnd_or <- I_cnd_or | X[,Who] == Or_Targets[i]
+ }
+ I_cnd <- I_cnd & I_cnd_or
}
}
-
- I_jnt=I_mrg & I_cnd
-
- if ( !isempty( View[[k]]$Cond_Who ) )
- {
- New_A = View[[k]]$sgn %*% t( (I_jnt - View[[k]]$v * I_cnd) )
- New_b = 0
+
+ I_jnt <- I_mrg & I_cnd
+
+ if (!isempty(View[[k]]$Cond_Who)) {
+ New_A <- View[[k]]$sgn %*% t((I_jnt - View[[k]]$v * I_cnd))
+ New_b <- 0
}
- else
- {
- New_A = View[[k]]$sgn %*% t( I_mrg )
- New_b = View[[k]]$sgn %*% View[[k]]$v
+ else {
+ New_A <- View[[k]]$sgn %*% t(I_mrg)
+ New_b <- View[[k]]$sgn %*% View[[k]]$v
}
-
- A = rbind( A , New_A ) # constraint for the conditional expectations...
- b = rbind( b , New_b )
- g = rbind( g , -log( 1 - View[[k]]$c ) )
+
+ A <- rbind(A, New_A) # constraint for the conditional expectations...
+ b <- rbind(b, New_b)
+ g <- rbind(g, -log(1 - View[[k]]$c))
}
- return( list( A = A , b = b , g = g ) )
+ return(list(A = A, b = b, g = g))
}
#' tweak a matrix
-#' @param A matrix A consisting of inequality constraints ( Ax <= b )
-#' @param b matrix b consisting of inequality constraint vector b ( Ax <= b )
-#' @param g TODO: TBD
+#' @param A matrix A consisting of inequality constraints (Ax <= b)
+#' @param b matrix b consisting of inequality constraint vector b (Ax <= b)
+#' @param g TODO: TBD
#'
#' @return db
#' @export
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
-Tweak = function( A , b , g )
-{
- library( matlab )
- library( limSolve )
-
- K = nrow( A )
- J = ncol( A )
-
- g_ = rbind( g , zeros( J , 1 ) )
-
- Aeq_ = cbind( zeros( 1 , K ) , ones( 1 , J ) )
- beq_ = 1
-
- lb_ = rbind( zeros( K , 1 ) , zeros( J , 1 ) )
- ub_ = rbind( Inf * ones( K , 1 ) , ones( J , 1 ) )
-
- A_ = cbind( -eye( K ) , A )
- b_ = b
-
+
+Tweak <- function(A, b, g) {
+ library(matlab)
+ library(limSolve)
+
+ K <- nrow(A)
+ J <- ncol(A)
+
+ g_ <- rbind(g, zeros(J, 1))
+
+ Aeq_ <- cbind(zeros(1, K), ones(1, J))
+ beq_ <- 1
+
+ lb_ <- rbind(zeros(K, 1), zeros(J, 1))
+ ub_ <- rbind(Inf * ones(K, 1), ones(J, 1))
+
+ A_ <- cbind(-eye(K), A)
+ b_ <- b
+
# add lower-bound and upper-bound constraints
- A_ = rbind( A_ , -eye(ncol(A_)) )
- b_ = rbind( b_ , zeros( ncol(A_), 1) )
-
- x0 = rep( 1/ncol( Aeq_ ) , ncol( Aeq_ ) )
- # db_ = linprog( g_ , A_ , b_ , Aeq_ ,beq_ , lb_ , ub_ ) # MATLAB version
- # optimResult = linp( E = Aeq_ , # matrix containing coefficients of equality constraints Ex=F
- # F = beq_ , # vector containing the right-hand side of equality constraints
- # G = -1*A_ , # matrix containint coefficients of the inequality constraints GX >= H
- # H = -1*b_ , # vector containing the right-hand side of the inequality constraints
- # Cost = -1*g_ , # vector containing the coefficients of the cost function
- # ispos = FALSE )
-
- costFunction = function( x ) { matrix( x , nrow = 1 ) %*% matrix( -1*g_ , ncol = 1) }
- gradient = function( x ) { -1*g_ }
- optimResult = optim( par = x0 ,
- fn = costFunction , # CHECK
- gr = gradient ,
+ A_ <- rbind(A_, -eye(ncol(A_)))
+ b_ <- rbind(b_, zeros(ncol(A_), 1))
+
+ x0 <- rep(1/ncol(Aeq_), ncol(Aeq_))
+
+ # db_ = linprog(g_, A_, b_, Aeq_,beq_, lb_, ub_) # MATLAB version
+ # matrix containing coefficients of equality constraints Ex=F
+ # optimResult = linp(E = Aeq_,
+ # vector containing the right-hand side of equality constraints
+ # F = beq_,
+ # matrix containint coefficients of the inequality constraints GX >= H
+ # G = -1*A_,
+ # vector containing the right-hand side of the inequality constraints
+ # H = -1*b_,
+ # vector containing the coefficients of the cost function
+ # Cost = -1*g_,
+ # ispos = FALSE)
+
+ costFunction <- function(x) {
+ matrix(x, nrow = 1) %*% matrix(-1*g_, ncol = 1)
+ }
+ gradient <- function(x) {
+ -1*g_
+ }
+ optimResult <- optim(par = x0,
+ fn = costFunction, # CHECK
+ gr = gradient,
method = "L-BFGS-B",
- lower = lb_ ,
- upper = ub_ ,
- hessian = FALSE )
-
- # library( linprog )
- # optimResult2 = solveLP( E = Aeq_ , # numeric matrix containing coefficients of equality constraints Ex=F
- # F = beq_ , # numeric vector containing the right-hand side of equality constraints
- # G = -1*A_ , # numeric matrix containint coefficients of the inequality constraints GX >= H
- # H = -1*b_ , # numeric vector containing the right-hand side of the inequality constraints
- # Cost = -g_ , # numeric vector containing the coefficients of the cost function
- # ispos = FALSE )
-
- db_ = optimResult$X
-
- db = db_[ 1:K ]
-
- return( db )
+ lower = lb_,
+ upper = ub_,
+ hessian = FALSE)
+
+ # library(linprog)
+ # optimResult2 = solveLP(E = Aeq_,
+ # vector containing the right-hand side of equality constraints
+ # F = beq_,
+ # matrix containint coefficients of the inequality constraints GX >= H
+ # G = -1*A_,
+ # vector containing the right-hand side of the inequality constraints
+ # H = -1*b_,
+ # vector containing the coefficients of the cost function
+ # Cost = -1*g_,
+ # ispos = FALSE)
+
+ db_ <- optimResult$X
+
+ db <- db_[1:K]
+
+ return(db)
}
-#' Takes a matrix of joint-scenario probability distributions and generates expectations, standard devation, and correlation matrix for the assets
+#' @title Takes a matrix of joint-scenario probability distributions and
+#' generates expectations, standard devation, and correlation matrix
+#' for the assets
#'
-#' Takes a matrix of joint-scenario probability distributions and generates expectations, standard devation, and correlation matrix for the assets
+#' @description Takes a matrix of joint-scenario probability distributions and
+#' generates expectations, standard devation, and correlation
+#' matrix for the assets
#'
-#' @param X a matrix of joint-probability scenarios (rows are scenarios, columns are assets)
-#' @param p a numeric vector containing the probabilities for each of the scenarios in the matrix X
+#' @param X a matrix of joint-probability scenarios (rows are scenarios,
+#' columns are assets)
+#' @param p a numeric vector containing the probabilities for each of the
+#' scenarios in the matrix X
#'
-#' @return means a numeric vector of the expectations (probability weighted) for each asset
-#' @return sd a numeric vector of standard deviations corresponding to the assets in the covariance matrix
-#' @return correlationMatrix the correlation matrix resulting from converting the covariance matrix to a correlation matrix
+#' @return means a numeric vector of the expectations
+#' (probability weighted) for each asset
+#' @return sd a numeric vector of standard deviations
+#' corresponding to the assets in the covariance
+#' matrix
+#' @return correlationMatrix the correlation matrix resulting from
+#' converting the covariance matrix to a
+#' correlation matrix
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
#' @export
-ComputeMoments = function( X , p )
-{
- library( matlab )
- J = nrow( X ) ; N = ncol( X )
- m = t(X) %*% p
- Sm = t(X) %*% (X * repmat( p , 1 , N ) ) # repmat : repeats/tiles a matrix
- S = Sm - m %*% t(m)
- C = cov2cor( S ) # the correlation matrix
- s = sqrt( diag( S ) ) # the vector of standard deviations
- return( list( means = m , sd = s , correlationMatrix = C ) )
+
+ComputeMoments <- function(X, p) {
+ library(matlab)
+ N <- ncol(X)
+ m <- t(X) %*% p
+ Sm <- t(X) %*% (X * repmat(p, 1, N)) # repmat : repeats/tiles a matrix
+ S <- Sm - m %*% t(m)
+ C <- cov2cor(S) # the correlation matrix
+ s <- sqrt(diag(S)) # the vector of standard deviations
+ return(list(means = m, sd = s, correlationMatrix = C))
}
Modified: pkg/Meucci/R/HermiteGrid.R
===================================================================
--- pkg/Meucci/R/HermiteGrid.R 2015-06-18 22:35:26 UTC (rev 3681)
+++ pkg/Meucci/R/HermiteGrid.R 2015-06-20 08:11:15 UTC (rev 3682)
@@ -1,67 +1,74 @@
#' Generates the normalized probability for an input probability value
#'
-#' @param p a numeric value containing the probability value required to be normalized
+#' @param p a numeric value containing the probability value
+#' required to be normalized
#'
-#' @return normalizedp a numeric value containing the normalized probability value
+#' @return normalizedp a numeric value containing the normalized probability
+#' value
#'
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
#' @export
-normalizeProb = function( p )
-{
- tol = 1e-20
- tmp = p
- tmp[ tmp < tol ] = tol
- normalizedp = exp( log( tmp ) - log( sum( tmp ) ) )
-
- return( normalizedp )
+
+normalizeProb <- function(p) {
+
+ tol <- 1e-20
+ tmp <- p
+ tmp[tmp < tol] <- tol
+ normalizedp <- exp(log(tmp) - log(sum(tmp)))
+
+ return(normalizedp)
}
#' Generate the intervals containing jth point of the grid.
#'
-#' @param x a vector containing the scenarios
+#' @param x a vector containing the scenarios
#'
#' @return a list containing
#' @return xLB a vector containing the lower bound of each interval
#' @return xUB a vector containing the upper bound of each interval
-#"
-#' \deqn{ I_{j} \equiv \big[ x_{j} - \frac{x_{j} - x_{j-1} }{2}, x_{j} + \frac{x_{j+1} - x_{j}}{2} \big) }
-#' @references
+#'
+#' \deqn{ I_{j} \equiv \big[x_{j} - \frac{x_{j} - x_{j-1} }{2}, x_{j} +
+#' \frac{x_{j+1} - x_{j}}{2} \big) }
+#' @references
#' A. Meucci - "Fully Flexible Extreme Views" - See formula (17)
#' \url{http://ssrn.com/abstract=1542083}
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
#' @export
-subIntervals = function( x )
-{
- n = nrow( x )
- x = t(x)
- xMesh = rep( NaN , n + 1)
- xMesh[ 1 ] = x[ 1 ]
- xMesh[ n + 1 ] = x[ n ]
- xMesh[ 2:n ] = x[2:n] - 0.5 * ( x[ 2:n ] - x[ 1:n-1 ] ) # simple product
-
- # cadlag mesh
- xUB = xMesh[ -1 ] - 2.2e-308 # right
- xLB = xMesh[ -(n + 1) ] # left
-
- return( list( xLB = xLB , xUB = xUB ) )
+
+subIntervals <- function(x) {
+
+ n <- nrow(x)
+ x <- t(x)
+ xMesh <- rep(NaN, n + 1)
+ xMesh[1] <- x[1]
+ xMesh[n + 1] <- x[n]
+ xMesh[2:n] <- x[2:n] - 0.5 * (x[2:n] - x[1:n - 1]) # simple product
+
+ # cadlag mesh
+ xUB <- xMesh[-1] - 2.2e-308 # right
+ xLB <- xMesh[-(n + 1)] # left
+
+ return(list(xLB = xLB, xUB = xUB))
}
-#' Integrate the subinterval for the given cumulative distribution function to get the equivalent probability
+#' Integrate the subinterval for the given cumulative distribution function to
+#' get the equivalent probability
#'
#' @param x a vector containing the data points
#' @param cdf the cumulative distribution function
#'
-#' @return p a vector containing the cdf evaluated for each of the subintervals
-#'
+#' @return p a vector containing the cdf evaluated for each of
+#' the subintervals
#' @author Ram Ahluwalia \email{ram@@wingedfootcapital.com}
#' @export
-integrateSubIntervals = function( x , cdf )
-{
- bounds = subIntervals( x )
-
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/returnanalytics -r 3682
More information about the Returnanalytics-commits
mailing list