1 files changed, 50 insertions, 29 deletions
diff --git a/R/distrib.R b/R/distrib.R
index 6c588ff..edc5081 100644
--- a/R/distrib.R
+++ b/R/distrib.R
@@ -15,8 +15,8 @@
 #' \code{GHquad} computes the quadrature weights for integrating against a
 #' Gaussian distribution.
 #'
-#' if f is a function, then with(GHquad(100), crossprod(f(Z), w))
-#' will compute \eqn{\frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty f(x)e^{-\frac{x^2}{2}}\,dx}.
+#' if f is a function, then \eqn{\sum_{i=1}^n f(Z_i)w_i \approx
+#' \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty f(x)e^{-\frac{x^2}{2}}\,dx}.
 #' @param n Integer, the number of nodes
 #' @return A list with two components:
 #'  \item{Z}{the list of nodes}
@@ -36,11 +36,11 @@ GHquad <- function(n){
 #' \code{lossdistrib} computes the probability distribution of a sum
 #' of independent Bernouilli variables with unequal probabilities.
 #'
-#' This uses the basic recursive algorithm of Andersen, Sidenius and Basu
-#' We compute the probability distribution of S = \sum_{i=1}^n X_i
-#' where X_i is Bernouilli(p_i)
+#' We compute the probability distribution of \eqn{S = \sum_{i=1}^n X_i}
+#' where \eqn{X_i} is Bernouilli(\eqn{p_i}). This uses the simple recursive
+#' algorithm of Andersen, Sidenius and Basu
 #' @param p Numeric vector, the vector of success probabilities
-#' @return A vector q such that q[k]=P(S=k)
+#' @return A vector q such that \eqn{q_k=\Pr(S=k)}
 lossdistrib <- function(p){
     ## basic recursive algorithm of Andersen, Sidenius and Basu
     n <- length(p)
@@ -58,13 +58,13 @@ lossdistrib <- function(p){
 #' \code{lossdistrib.fft} computes the probability distribution of a sum
 #' of independent Bernouilli variables with unequal probabilities.
 #'
-#' This uses the fft. Complexity is of order O(n m) + O(m\log{m})
-#' where m is the size of the grid and n, the number of probabilities.
+#' We compute the probability distribution of \eqn{S = \sum_{i=1}^n X_i}
+#' where \eqn{X_i} is Bernouilli(\eqn{p_i}).
+#' This uses the FFT, thus omplexity is of order \eqn{O(n m) + O(m\log(m))}
+#' where \eqn{m} is the size of the grid and \eqn{n}, the number of probabilities.
 #' It is slower than the recursive algorithm in practice.
-#' We compute the probability distribution of S = \sum_{i=1}^n X_i
-#' where X_i is Bernouilli(p_i)
 #' @param p Numeric vector, the vector of success probabilities
-#' @return A vector such that q[k]=P(S=k)
+#' @return A vector such that \eqn{q_k=\Pr(S=k)}
 lossdistrib.fft <- function(p){
     n <- length(p)
     theta <- 2*pi*1i*(0:n)/(n+1)
@@ -73,13 +73,20 @@ lossdistrib.fft <- function(p){
     return(1/(n+1)*Re(fft(v)))
 }
 
-#' recursive algorithm with first order correction
+#' Loss distribution of a portfolio
+#'
+#' \code{lossdistrib2} computes the probability distribution of a sum
+#' of independent Bernouilli variables with unequal probabilities.
 #'
+#' We compute the probability distribution of \eqn{L = \sum_{i=1}^n w_i S_i X_i}
+#' where \eqn{X_i} is Bernouilli(\eqn{p_i}). If \code{defaultflag} is TRUE, we
+#' compute the distribution of \eqn{D = \sum_{i=1}^n w_i X_i} instead.
+#' This a recursive algorithm with first order correction for discretization.
 #' @param p Numeric, vector of default probabilities
 #' @param w Numeric, vector of weights
 #' @param S Numeric, vector of severities
 #' @param N Integer, number of ticks in the grid
-#' @param defaultflag Boolean, if True, we compute the default distribution
+#' @param defaultflag Boolean, if TRUE, we compute the default distribution
 #' (instead of the loss distribution).
 #' @return a Numeric vector of size \code{N} computing the loss (resp.
 #' default) distribution if \code{defaultflag} is FALSE (resp. TRUE).
@@ -106,12 +113,15 @@ lossdistrib2 <- function(p, w, S, N, defaultflag=FALSE){
     return(q)
 }
 
-#' recursive algorithm with first order correction truncated version
-#' this is actually slower than lossdistrib2. But in C this is
-#' twice as fast.
-#' For high severities, M can become bigger than N, and there is
-#' some probability mass escaping.
+#' Loss distribution truncated version
+#'
+#' \code{lossdistrib2.truncated} computes the probability distribution of a sum
+#' of independent Bernouilli variables with unequal probabilities up
+#' to a cutoff N.
 #'
+#' This is actually slower than \code{lossdistrib2}, but in C this is
+#' twice as fast. For high severities, M can become bigger than the cutoff, and
+#' there is some probability mass escaping.
 #' @param p Numeric, vector of default probabilities
 #' @param w Numeric, vector of weights
 #' @param S Numeric, vector of severities
@@ -140,15 +150,26 @@ lossdistrib2.truncated <- function(p, w, S, N, cutoff=N){
     return(q)
 }
 
+#' Recovery distribution of a portfolio
+#'
+#' \code{recovdist} computes the recovery distribution of portfolio
+#' described by a vector of default probabilities, and prepay probabilities.
+#' \eqn{R=\sum_{i=1}^n w_i X_i} where \eqn{X_i=0} w.p. \eqn{1-dp_i-pp_i},
+#' \eqn{X_i=1-S_i} with probability \eqn{dp_i}, and \eqn{X_i=1} w.p. \eqn{pp_i}
+#'
+#' It is a recursive algorithm with first-order correction. For a unit of loss
+#' \eqn{lu}, each non-zero value \eqn{v} is interpolated on the grid
+#' as the pair of values
+#' \eqn{\left\lfloor\frac{v}{lu}\right\rfloor}  and
+#' \eqn{\left\lceil\frac{v}{lu}\right\rceil} so that \eqn{X_i} has
+#' four non zero values.
+#' @param dp Numeric, vector of default probabilities
+#' @param pp Numeric, vector of prepay probabilities
+#' @param w Numeric, vector of weights
+#' @param S Numeric, vector of severities
+#' @param N Integer, number of ticks in the grid
+#' @return a Numeric vector of size \code{N} computing the recovery distribution
 recovdist <- function(dp, pp, w, S, N){
-    ## computes the recovery distribution for a sum of independent variables
-    ## R=\sum_{i=1}^n w[i] X_i
-    ## where X_i = 0             w.p 1 - dp[i] - pp[i]
-    ##           = 1 - S[i]      w.p dp[i]
-    ##           = 1             w.p pp[i]
-    ## each non zero value v is interpolated on the grid as
-    ## the pair of values floor(v/lu) and ceiling(v/lu) so that
-    ## X_i has four non zero values
     n <- length(dp)
     q <- rep(0, N)
     q[1] <- 1
@@ -173,16 +194,16 @@ recovdist <- function(dp, pp, w, S, N){
     return(q)
 }
 
-#' recursive algorithm with first order correction to compute the joint
+#' Joint loss-recovery distributionrecursive algorithm with first order correction to compute the joint
 #' probability distribution of the loss and recovery.
+#'
 #' For high severities, M can become bigger than N, and there is
 #' some probability mass escaping.
-#'
 #' @param p Numeric, vector of default probabilities
 #' @param w Numeric, vector of weights
 #' @param S Numeric, vector of severities
 #' @param N Integer, number of ticks in the grid
-#' @param cutoff Integer, where to stop computing the exact probabilities
+#' @param defaultflab Logical, whether to return the loss or default distribution
 #' @return q Matrix of joint loss, recovery probability distribution
 #' colSums(q) is the recovery distribution marginal
 #' rowSums(q) is the loss distribution marginal