3
$\begingroup$

I want to compute the variance of $FV =A\sum_{k=1}^n(1+r)^k=:g(r)$, assuming interest rate $r\sim\mathrm{N}(\mu,\sigma^2)$ and constant equal payments $A$. As discussed in the comments, the Delta Method gives an approximation $\mathrm{Var}[FV]≈\sigma^2\times(g′(\mu))^2$. Is there a way to get an exact answer?

$\endgroup$
3
  • $\begingroup$ Why is the (random) interest the same in every period? Do you really mean for the daily return - not it's distribution, but its actual realized value - to be exactly the same for day 1 and day 2 and day 3? If not, you need a different symbol for each day ($r_1,r_2,r_3,...$) $\endgroup$
    – Glen_b
    Commented Jun 6, 2015 at 23:39
  • $\begingroup$ Thank you @Zen for the formula. I compared the distribution with $$ Var(FV)≈\sigma^2×(g′(μ))2. $$ and with a monte carlo simulation. The spread is way smaller with the approximation, so I guess there's a tradeoff between accuracy and speed. $\endgroup$
    – CCL
    Commented Jun 7, 2015 at 11:46
  • 1
    $\begingroup$ Way better !! Inputs : pension of € 3000/month during 240 months, starting at now+40 years Output (P(X<=x)=20%, P(X<=x)=50%, P(X>=x)=20%) : Formula : € 2,506,743.83 € 2,513,886.71 € 2,521,029.58 MonteCarlo : € 2,506,550.77 € 2,519,574.13 € 2,532,597.50 Thanks, I'm going to edit my answer with your formula $\endgroup$
    – CCL
    Commented Jun 8, 2015 at 13:52

1 Answer 1

4
$\begingroup$

Let $X = 1+r$ have a Normal$(\mu,\sigma)$ distribution (so that the $\mu$ of this answer is $1$ more than the $\mu$ of the question). The question concerns

$$\eqalign{ \text{Var}(Y) &= A^2 \text{Var}\left(\sum_{i=1}^n X^i\right) = A^2\sum_{i,j=1}^n \text{Cov}(X^i,X^j) \\ &= A^2\sum_{i,j=1}^n \left(\mathbb{E}(X^{i+j}) - \mathbb{E}(X^i)\mathbb{E}(X^j)\right) }$$

where $Y = \sum_{i=1}^n A X^i$. This can be computed as a simple, fast outer product based on the vector $\nu = (\nu_1,\nu_2, \ldots, \nu_{2n})$ of expectations

$$\nu_i = \mathbb{E}(X^i) = \begin{cases} \sigma^k(k-1)\text{!!} \, _1F_1\left(-\frac{k}{2};\frac{1}{2};-\frac{\mu ^2}{2\sigma^2}\right) & k \text{ even} \\ \sigma^k\frac{\mu}{\sigma}\, k\text{!!} \, _1F_1\left(\frac{1-k}{2};\frac{3}{2};-\frac{\mu ^2}{2\sigma^2}\right) & \text{ odd} \end{cases}.$$

Here, $i!! = i(i-2)(i-4)\cdots(3)(1)$ for any odd positive integer $i$ and $_1F_1$ is a polynomial equal to the confluent hypergeometric function

$$_1F_1(-j;(2b+1)/2;-z/2) = \sum_{i=0}^j \frac{(-j)_{(i)} (-z/2)^i}{i! ((2b+1)/2)_{(i)}}$$

with $(x)_{(i)} = (x)(x+1)\cdots(x+i-1)$.

This allows for computing the variance exactly in $O(n^2)$ basic arithmetic operations.


Code

The following R code implements this algorithm as the function v (which performs the outer product) and tests it with a realistic simulation (240 months at 3% annual rate with a monthly standard deviation of 0.2% in that rate). Specifically, v(n, 1+mu, sigma) is the variance of the sum of returns after n periods when the per-period mean of $r$ is mu and the per-period standard deviation of $r$ is sigma. The simulation agrees with the formula to within a fraction of a standard error.

(The caller is responsible for multiplying the output of v by $A^2$.)

Note that care must be taken to avoid overflow: extremely large values of the hypergeometric function will balance extremely small values of $\sigma^k$. The functions therefore work with logarithms throughout.

hypergeom.confluent.log <- function(a,b,z,offset) {
  # Special values of the log of the confluent hypergeometric function 1F1:
  # a is a negative integer; b is 1/2 or 3/2; z is negative.
  i <- 0:(-a)
  if (b == 1/2) {
    u <- offset + lfactorial(i) - lfactorial(2*i)
  } else {u <- log(2) + offset + lfactorial(i+1) - lfactorial(2*(i+1))}

  y <- lchoose(-a, i) + i*(log(-z) + log(4)) + u
  y.max <- max(y)
  log(sum(exp(y - y.max))) + y.max
}

e <- function(k, mu, sigma) { 
  # Expectation of X^k where X has a Normal(mu, sigma) distribution.
  ff <- function(i) lfactorial(i) - lfactorial((i-1)/2) - (i-1)/2 * log(2)

  exp(k * log(sigma) + ifelse (k %% 2 == 0, 
      hypergeom.confluent.log(-k/2, 1/2, -mu^2/(2*sigma^2), ff(k-1)),
      log(mu / sigma) + hypergeom.confluent.log((1-k)/2, 3/2,-mu^2/(2*sigma^2), ff(k))))
}

v <- function(n, mu, sigma) {
  # Variance of X + X^2 + ... + X^n where X has a Normal(mu, sigma) distribution.
  x <- sapply(1:(2*n), function(k) e(k, mu, sigma))
  sum(matrix(x[outer(1:n, 1:n, function(i,j) i+j)], n) - outer(x[1:n], x[1:n]))
}

f <- function(n, r) {
  # Sum of (1+r) + (1+r)^2 + ... + (1+r)^n
  x <- 1+r
  (x^(n+1) - x) / r
}
#
# Simulation test of `v`.
#
n <- 300       # Periods
mu <- .03/12   # Expected interest rate per period
sigma <- 0.002 # SD of periodic interest rate
n.sim <- 1e6   # Simulation size
set.seed(17)
sim <- sapply(rnorm(n.sim, mu, sigma), function(r) f(n, r))
hist(sim)
#
# Compare the variance of the simulated results to the formula.
#
m <- mean(sim)
m.4 <- mean((sim-m)^4)
m.2 <- mean((sim-m)^2)
se <- sqrt(m.4/n.sim - m.2^2*(n.sim-3) / (n.sim*(n.sim-1)))
v.ref <- v(n, 1+mu, sigma)
v.sim <- var(sim)
c(Variance=v.sim, SE=se, Formula=v.ref, Z=(v.sim - v.ref) / se)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.