Archive for the ‘R’ Category

Mountain Rainfall

Minor shifts in atmospheric circulation can have major impacts on rainfall patterns in mountainous regions. This has important implications for crops such as coffee which are grown at altitude.

Rainfall derives from uplift of moist air. In mountainous areas, uplift is forced by terrain and so depends on the horizontal wind vector \mathbf U. The dependence of rainfall on wind direction is illustrated below for a fictitious range of hills (dashed contours indicate the hills).

The “upslope” model for rainfall (P) at location \mathbf{r} is:

(1)   \begin{equation*} P(\mathbf{r}) = P_o + \rho_s  \mathbf{U} \cdot  \vec \nabla h\left(\mathbf{r}\right) \end{equation*}

\vec\nabla h is the gradient of the terrain height and \rho_s is related to the saturated water vapour density. P_0 is the background rainfall rate.

In reality, Equation (2) is a poor approximation.  Rain spreads downwind because there is a delay before raindrops can form, and another delay before fallout. Airflow dynamics also spreads the influence of an obstacle upwind. These effects can be included in a Fourier transformed version of equation (1):

(2)   \begin{equation*} P(\mathbf{r}) = P_0 + \rho_s \sum_\mathbf{k} g(\mathbf{k}) exp\left( i \mathbf{k}\cdot \mathbf{r}\right) i \mathbf{U} \cdot \mathbf{k} \tilde h(\mathbf{k}) \end{equation*}

\tilde h(\mathbf{k}) is the Fourier transform of the terrain height and g(\mathbf{k}) is a form factor which includes the effects of advection and mountain wave dynamics. Equation (2) reduces to the upslope model when g(\mathbf{k})=1.

R: Fast Fourier Transform & Rcpp

The figure illustrates orographic rainfall for a fictitious range of gaussian hills. It is very quick to compute (2) using Discrete FFT. If you use R, things can be speeded up by constructing the integrand in C++ using Rcpp before passing it to fft. This is very fast even for realistic landscapes.

C++/Rcpp implementation


Hidden history of the bond market

Sometimes time-series data derive from an underlying system which can exist in multiple distinct states. Hidden Markov Models (HMMs) are a popular approach in this situation. In a particular state, the data fluctuate with variance characteristic of that state. An N x N transition probability matrix describes the additional dynamics associated with flipping between the N states. The value of N and nature of the unobserved states may or may not be obvious depending on the problem.

Yields on long-term British government debt for the period 1727-2013 are available from the UK debt management office (annual mean yields on perpetual bonds or “consols”). Bond yields change in response to market perceptions of risk. Historically, volatile periods tend to be associated with large attritional wars, such as Seven Years War (1754-1763), Napoleonic wars (1803-1815), World War One (1914-1918) etc.

It is tempting to fit a HMM to the (differenced) gilt yield data. N is found by minimising the Bayesian Information Criterion (BIC). It turns out that N=3 is optimal for the differenced time-series (code and graph below).

The modern era of high inflation requires a third underlying state “I” which lasted from 1968-1999.


The R code below returns N=3.

library(RHmm) <- “”
gilt <- read.csv(
returns <- diff(gilt$yield)
bics <- sapply(2:10, function(nn) HMMFit(returns, dis=”NORMAL”, nStates=nn)$BIC)