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?
-
$\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_bCommented 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$– CCLCommented 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$– CCLCommented Jun 8, 2015 at 13:52
1 Answer
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)