Variance of the future value of a cash flow 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?
 A: 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)

