{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Markdown parser included in IPython is MathJax-aware. This means that you can freely mix in mathematical expressions using the MathJax subset of Tex and LaTeX"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical properties of the rainfall model\n",
    "-------------------------------------------------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In model fitting it is usually necessary to use equations for aggregated properties, because rainfall data are usually sampled over discrete time intervals. Let $Y_{ij}^{h}(x)$ be the aggregated time series of rainfall due to type $i$ storms at point $x = (x, y) ∈ R^2$ in the $j-th$ time interval of width $h$, and let $Y_{j}^h(x)$ be the total rainfall in the $j-th$ interval due to the superposition of the $n$ storm types [1]. Then,\n",
    "\n",
    "$$Y_{j}^h (x) = \\sum_{i=1}^{n} \\int_{(j-1)/h}^{jh} Y_i (x, t)dt$$\n",
    "\n",
    "where $Y_i (x, t)$ is the rainfall intensity at point $x$ and time $t$ due to type $i$ storms $(i = 1, ..., n)$. Since the superposed processes are independent, statistical properties of the aggregated time series follow just by summing the various properties given below [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Mean ([2] Eq.12):**\n",
    "\n",
    "$$\\mu(h)= \\sum_{i=1}^{n} E\\lbrace Y_{ij}^{h}(x)\\rbrace=h\\sum_{i=1}^{n}\\frac{\\lambda_{i}\\nu_{i}}{\\chi_{i}\\eta_{i}}$$\n",
    "\n",
    "in the STNSRP model, $\\chi_{i}$ must me change by $\\chi_{i}^-1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Covariance ([2] Eq.14):**\n",
    "\n",
    "$$\\gamma(x,x, l, h) =  \\sum_{i=1}^{n} Cov\\lbrace Y_{ij}^{h}, Y_{i,j+l}^{h}\\rbrace = \\frac{\\lambda_i(\\nu_{i}^{2}-1)[\\beta_i^{3}A_i(h, l)-\\eta_i^{3}B_i(h,l)]}{\\beta_i\\xi^{2}\\eta_i^{3}(\\beta_i^{2}-\\eta_i^{2})} - \\frac{4\\lambda_i \\chi_{i}A_i(h,l)}{\\xi_i^{2}\\eta_i^{3}}$$\n",
    "\n",
    "when $l =0$,\n",
    "\n",
    "$$A_i(h,l)=A_i(h)=\\eta_i h -1 + e^{-\\eta_i h}$$\n",
    "$$B_i(h,l)=B_i(h)=\\beta_i h -1 + e^{-\\beta_i h}$$\n",
    "\n",
    "when $l >0$,\n",
    "\n",
    "$$A_i(h,l)=0.5(1-e^{-\\eta_i h})^{2}e^{-\\eta_i h(l-1)}$$\n",
    "$$B_i(h,l)=0.5(1-e^{-\\beta_i h})^{2}e^{-\\beta_i h(l-1)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Probability of not rain in an arbitrary time of length $h$ ([3] Eq.6):**\n",
    "\n",
    "The probability that an arbitrary time interval $[(j − 1)h, j h]$ is dry at a point is obtained by multiplying\n",
    "the probabilities of the independent processes and is given by the following:\n",
    "    \n",
    "$$\\phi(h) = \\sum_{i=1}^{n} Pr\\lbrace Y_{ij}^{(h)}(x)=0 \\rbrace = exp(-\\lambda_i h + \\lambda_i \\beta_i^{-1}(\\nu_{i} -1)^{-1}\\lbrace1-exp[1-\\nu_{i} + (\\nu_{i}-1) e^{-\\beta_i h}]\\rbrace -\\lambda_i\\int_{0}^{\\infty}[1-p_{h}(t)]dt)$$\n",
    "    \n",
    "    \n",
    "We use the approximation shown in [2] Eq.17 to avoid having to solve the integral:\n",
    "\n",
    "$$\\int_{0}^{\\infty}[1-p_{h}(t)]dt=\\beta_i^{-1}[\\gamma +ln(\\alpha \\nu_{i}-1)]$$\n",
    "\n",
    "where \\\\(\\gamma =0.5771 \\\\) y \\\\(\\alpha_i = \\eta_i/(\\eta_i - \\beta_i)-e^{-\\beta_i h}\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Transition probabilities**:\n",
    "\n",
    "The transition probabilities, $pr \\lbrace Y_{i,j+1}^{(h)}(x)>0 | Y_{ij}^{(h)}(x)>0 \\rbrace$ and $pr \\lbrace Y_{i,j+1}^{(h)}(x)=0 | Y_{i}^{(h)}(x)=0 \\rbrace$, denoted as $\\phi_{WW}(h)$ y $\\phi_{DD}(h)$, respectively, can be expressed in terms of \\\\(\\phi(h)\\\\) following ([3] Eq.7,8 and 9):\n",
    "\n",
    "$$\\phi_{DD}(h)=\\phi(2h)/\\phi(h)$$\n",
    "\n",
    "$$\\phi(h)=\\phi_{DD}(h) + \\lbrace 1-\\phi_{WW}(h) \\rbrace \\lbrace(1-\\phi(h))\\rbrace$$\n",
    "\n",
    "$$\\phi_{WW}(h)=\\lbrace 1-2 \\phi(h) + \\phi(2h) \\rbrace\\lbrace 1 -\\phi(h) \\rbrace$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The third moment function ([4] Eq.10):**\n",
    "\n",
    "$$\\xi_{h}=E\\lbrace Y_{j}^{(h)}(x)-\\mu(h)\\rbrace^{3}= \\sum_{i=1}^{n}[ 6\\lambda_i\\nu_{i}E(X^{3})(\\eta_i h -2 +\\eta_i h e^{-\\eta_i h}+2e^{-\\eta_i h})/\\eta_i^{4} + 3\\lambda_i\\chi_i E(X^{2})\\nu_{i}^2 f(\\eta_i, \\beta_i, h)/{\\lbrace2\\eta_i^{4}\\beta_i (\\beta_i^{2}-\\eta_i^{2})^{2}}\\rbrace + \\lambda_i \\chi_i^{3}\\nu_{i}^{3} g(\\eta_i, \\beta_i, h)/{\\lbrace e\\eta_i^{4}\\beta_i(e\\eta_i^{2}-\\beta_i^{2})(\\eta_i-\\beta_i)(2\\beta_i +\\eta_i)(\\beta_i +2\\eta_i)}\\rbrace]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the STNSRP model, C is a Poisson random variable, so that $E\\lbrace C(C-1)\\rbrace = \\nu^2$ and $E\\lbrace C(C-1)(C-2)\\rbrace = \\nu^3$\n",
    "If is a geometric then $E\\lbrace C(C-1)\\rbrace = 2\\nu^2(\\nu-1)$ and $E\\lbrace C(C-1)(C-2)\\rbrace = 6\\nu^2(\\nu-1)^{2}$\n",
    "\n",
    "For exponential cell intensities, $E(X_{ijk})$ and $E(X_{ijk}^{k})$ are replaced by $2\\chi_i^{2}$ and $6\\chi_i^{3}$\n",
    " respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$f(\\eta_i, \\beta_i, h)$ and $g(\\eta_i, \\beta_i, h)$ are derived bellow:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$f(\\eta_i, \\beta_i, h) = -2\\eta^{3}\\beta^{2}\\exp(-\\eta h)-2\\eta^{3}\\eta^{2}\\exp(-\\beta h)+\\eta^{2}\\beta^{3}\\exp(-2\\eta h)+2\\eta^{4}\\beta \\exp(-\\eta h)+2\\eta^{4}\\beta \\exp(-\\beta h)+2\\eta^{3}\\beta^{2}\\exp(-(\\eta+\\beta) h)-2\\eta^{4}\\beta\\exp(-(\\eta+\\beta)h)-8\\eta^{3}\\beta^{3} h+11 \\eta^{2}\\beta^{3}-2\\eta^{4}\\beta+2\\eta^{3}\\eta^{2}+4\\eta\\beta^{5} h+4\\eta^{5}\\beta h-7 \\beta^{5}-4\\eta^{5}+8\\beta^{5} \\exp(-\\eta h)-\\beta^{5}\\exp(-2\\eta h)\\\n",
    "    -2 h \\eta^{3}\\beta^{3}\\exp(-\\eta h)-12\\eta^{2}\\beta^{3}\\exp(-\\eta h)+2 h \\eta \\beta^{5}\\exp(-\\eta h)+4 \\eta^{5}\\exp(-\\beta h)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$g(\\eta_i, \\beta_i,h)=12\\eta^{5}\\beta\\exp(-\\beta h)+9\\eta^{4}\\beta^{2}+12\\eta \\beta^{5}\\exp(-\\eta h)+9 \\eta^{2}\\beta^{4}+12\\eta^{3}\\beta^{3}\\exp(-(\\eta+\\beta)h)-\\eta^{2}\\beta^{4}\\exp(-2\\eta h)-12 \\eta^{3}\\beta^{3}\\exp(-\\beta h)-9\\eta^{5}\\beta-9\\eta \\beta^{5}-3\\eta\\beta^{5}\\exp(-2\\eta h)-\\eta^{4}\\beta^{2}\\exp(-2\\beta h)-12\\eta^{3}\\beta^{3}\\exp(-\\eta h)+6\\eta^{5}\\beta^{2} h-10\\beta^{4}\\eta^{3}h+6\\beta^{5}\\eta^{2}h-10\\beta^{3}\\eta^{4} h+4\\beta^{6}\\eta h-8\\beta^{2}\\eta^{4}\\exp(-\\beta h)+4\\beta\\eta^{6} h+12\\beta^{3}\\eta^{3}-8\\beta^{4}\\eta^{2}\\exp(-\\eta h)-6\\eta^{6}-6\\beta^{6}-2\\eta^{6}\\exp(-2\\beta h)-2\\beta^{6}\\exp(-2\\eta h)+8\\eta^{6}\\exp(-\\beta h)+8\\beta^{6}\\exp(-\\eta h)-3\\beta\\eta^{5}\\exp(-2\\beta h)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The following aggregated statistics only affect to the STNSRM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each storm, the number of cells $\\nu$ that overlap a point in $R^2$ is a Poisson random variable with mean ([5] Eq.3):\n",
    "    \n",
    "$$\\nu_C = \\nu \\phi_c^2/(2\\pi) $$\n",
    "\n",
    "where $\\nu_C$ denotes the two-dimensional Poisson process (cells per $km^2$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The probability that a cell overlaps a point $x$ given that it overlapped a point $y$ a distance $d$ from $x$ ([5], Eq.9) \n",
    "\n",
    "\n",
    "$$ P(\\phi, d)\\approx \\frac{1}{30} \\sum_{i=1}^{4} \\lbrace 2f(\\frac{2 \\pi i}{20}) + 4f(\\frac{2\\pi i + \\pi}{20})\\rbrace - \\frac{1}{30} f(0) $$\n",
    "\n",
    "where \n",
    "\n",
    "$$f(y)=(\\frac{\\phi d}{2 cosy}+1)exp(\\frac{-\\phi d}{2cosy}), 0 \\le y < \\pi/2, f(\\pi/2)$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cross-correlation ([5], Eq.6):**\n",
    "    \n",
    "$$\\gamma(x,y,l,h)=\\sum_{i=1}^{n} Cov\\lbrace Y_{ij}^{h}(x), Y_{i,j+l}^{h}(y)\\rbrace  = \\sum_{i=1}^{n} [\\gamma_i(x,x,l,h)-2\\lambda_i\\lbrace1-P(\\phi_i, d)\\rbrace\\nu_{Ci}E(X^2)A_i(h,l)/\\eta_i^3]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1]: Cowpertwait, P., Ocio, D., Collazos, G., De Cos, O., Stocker, C. Regionalised spatiotemporal rainfall and temperature models for flood studies in the Basque Country, Spain (2013) Hydrology and Earth System Sciences, 17 (2), pp. 479-494.\n",
    "\n",
    "[2]: Cowpertwait, P.S.P. Further developments of the neyman‐scott clustered point process for modeling rainfall (1991) Water Resources Research, 27 (7), pp. 1431-1438.\n",
    "\n",
    "[3]: Cowpertwait, P.S.P., O'Connell, P.E., Metcalfe, A.V., Mawdsley, J.A. Stochastic point process modelling of rainfall. I. Single-site fitting and validation (1996) Journal of Hydrology, 175 (1-4), pp. 17-46. \n",
    "\n",
    "[4]: Cowpertwait, P.S.P. A Poisson-cluster model of rainfall: High-order moments and extreme values (1998) Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 454 (1971), pp. 885-898.\n",
    "\n",
    "\n",
    "[5]: Cowpertwait, P.S.P., Kilsby, C.G., O'Connell, P.E. A space-time Neyman-Scott model of rainfall: Empirical analysis of extremes (2002) Water Resources Research, 38 (8), pp. 6-1-6-14.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}