{ "cells": [ { "cell_type": "markdown", "id": "329b70ac-fb96-46f9-93a5-a935043d835d", "metadata": {}, "source": [ "
\n", " Climate of the Earth system\n", "

Prof. Dr. Markus Meier
\n", " Leibniz Institute for Baltic Sea Research Warnemünde (IOW)
\n", " E-Mail: markus.meier@io-warnemuende.de

\n", "
" ] }, { "cell_type": "markdown", "id": "519b554f-e050-4dc0-9cf9-c8dd075a8033", "metadata": {}, "source": [ "# Spectrum\n", "\n", "- = the Fourier transform of the auto-covariance function of the time series (presents the variance per frequency)\n", "- math. definition: Let $\\mathbf{X_t}$ be an ergodic weakly stationary stochastic process with auto-covariance function $\\gamma(\\tau), ~\\tau=0,\\pm1,...~$. Then the spectrum (or power spectrum) $\\Gamma$ of $\\mathbf{X_t}$ is the Fourier transform $\\cal{F}$ of the auto-covariance function $\\gamma$:\n", "```{math}\n", ":label: spectrum\n", "\\begin{align*}\n", "\\Gamma(\\omega) &= \\cal{F}\\{\\gamma\\}(\\omega)\\\\\n", "&= \\sum_{\\tau=-\\infty}^{\\infty} \\gamma(\\tau)e^{-2\\pi i\\tau\\omega},~~~~~\\forall \\omega \\in \\left[-\\frac{1}{2},~\\frac{1}{2} \\right]\n", "\\end{align*}\n", "```\n", "- note that the largest frequency that a time series with time step of 1.0 can resolve is $\\omega=\\frac{1}{2}$\n", "- using Euler's formular\n", "```{math}\n", "e^{i\\phi} = cos(\\phi) + i~sin(\\phi)\n", "```\n", "- we can rewrite [](#spectrum) to\n", "```{math}\n", "\\Gamma(\\omega) = \\gamma(0) + 2\\sum_{\\tau=1}^{\\infty}\\gamma(\\tau)cos(2\\pi\\tau\\omega)\n", "```\n", "- the spectrum can be interpreted as the covariance between the auto-correlation function and the cosine function at different frequencies\n", "- characteristics of the spectrum:\n", " 1. spectrum is continuous and differentiable in $[-\\frac{1}{2},~\\frac{1}{2}]$ (unlike the discrete time series or auto-correlation function)\n", " 2. spectrum describes the distribution of variance across time scales: $Var(\\mathbf{X_t}) = \\gamma(0) = 2\\int_0^{\\Sigma}\\Gamma(\\omega)~dw$, $\\Gamma(\\omega)$ is a variance per frequency\n", " 3. $\\gamma(\\tau) = \\int_{-\\frac{1}{2}}^{\\frac{1}{2}}\\Gamma(\\omega)e^{2\\pi i\\tau\\omega}~dw$\n", " 4. $\\left. \\frac{d}{d\\omega}\\Gamma(\\omega)\\right|_{\\omega=0} =0$, spectrum must be flat at long time scales for stationary processes\n", " 5. $\\Gamma(\\omega) \\propto \\mathbf{X}^2~~\\Rightarrow~~\\sigma(\\Gamma(\\omega))\\propto E(\\Gamma(\\omega))$, statistical uncertainty is proportional to expectation value (large for peaks), $E(\\mathbf{X})=k$, $Var(\\mathbf{X})=2k$ with k as the number of degrees of freedom\n", "\n", "```{figure} figures/L16/L16_1_SST.PNG\n", "---\n", "width: 50%\n", "---\n", "

Figure 1: Spectral estimates of SST time series in the northern North Atlantic. Left is with seasonal cycles and right is also with the mean seasonal cycle removed.

\n", "```\n", "```{figure} figures/L16/L16_2_spectrum_AR.PNG\n", "---\n", "width: 50%\n", "---\n", "

Figure 2: The spectrum of an AR(1)-process time series presented in log-linear scaling (upper) and in log-log scaling (lower). In addition the 10% and 90% quantilesof the spectral coefficients estimate are plotted (right).

\n", "```\n", "- for spectral variance per frequency on a log-linear scale as in the top row of Fig. 2, the area is the total variance. drawback: theoretical models follow power laws\n", "\n", "```{figure} figures/L16/L16_3_spectrum_elnino.PNG\n", "---\n", "width: 50%\n", "---\n", "

Figure 3: The spectrum of the monthly mean El Nino time series in four different presentations of the same spectrum. Note the lower right panel shows $\\omega \\Gamma(\\omega)$. For all spectra the fitted AR(1)-process spectrum with the 90% confidence interval is plotted for comparison.

\n", "```\n", "- peak at frequency of about $0.25yr^{-1}$, note the different representations of that peak.\n", "\n", "```{figure} figures/L16/L16_4_spectrum_ts.PNG\n", "---\n", "width: 50%\n", "---\n", "

Figure 4: Comparison of the time series (upper right), band-pass filtered time series (lower right) and its relation to the spectrum of the time series. The colored lines in the spectrum (left) illustrate the frequency-band for which the time series on the right were filtered. The numbers next to it indicate the mean spectral variance for the frequency-band.

\n", "```" ] }, { "cell_type": "markdown", "id": "1a0b40d6-0f08-475e-a57d-85f5715eeda8", "metadata": {}, "source": [ "## Spectra of AR(p)-processes\n", " - the spectrum of an AR(p)-process with process parameters $\\{\\alpha_1,...,\\alpha_p \\}$ and noise variance $Var(\\mathbf{Z_t})= \\sigma_Z^2$ is:\n", "```{math}\n", ":label: arp\n", "\\Gamma(\\omega) = \\frac{\\sigma_Z^2}{\\left| 1-\\sum_{k=1}^p \\alpha_ke^{-2\\pi ik\\omega} \\right|^2}\n", "```\n", "- spectrum of a white noise process AR(0): variances are equally distributed on all frequencies:\n", "```{math}\n", "\\Gamma(\\omega) = \\sigma_Z^2\n", "```\n", "#### Spectrum of an AR(1)-process:\n", "- for $p=1$, [](#arp) can be rewritten as:\n", "```{math}\n", ":label: spec\n", "\\Gamma(\\omega) = \\frac{\\sigma_Z^2}{\\left| 1-\\alpha_1e^{-2\\pi i\\omega} \\right|^2} = \\frac{\\sigma_Z^2}{1 + \\alpha_1^2-2\\alpha_1cos(2\\pi\\omega)}\n", "```\n", "- for small values $\\omega \\in [0,~\\frac{1}{2}]$ we can use the Taylor-expansion\n", "```{math}\n", ":label: exp\n", "cos(2\\pi\\omega) = 1- \\frac{(2\\pi\\omega)^2}{2!} + ...\n", "```\n", "- inserting the expansion [](#exp) into [](#spec) gives us the expression:\n", "```{math}\n", "\\Gamma(\\omega) \\approx \\frac{\\sigma_Z^2}{1 + \\alpha_1^2-2\\alpha_1(1-\\frac{(2\\pi\\omega)^2}{2!})} = \\frac{\\sigma_Z^2}{1 + \\alpha_1^2-2\\alpha_1 +\\alpha_1(2\\pi\\omega)^2} = \\frac{c_1\\sigma_Z^2}{c_2+\\omega^2}\n", "```\n", "- we can see that the spectrum follows a linear gradient of $-2$ in log-log scale for $\\omega >>0$\n", "- for $\\omega \\in ]0,~\\frac{1}{2}[$ the spectrum has no extrema because:\n", "```{math}\n", "\\frac{d}{d\\omega}\\Gamma(\\omega) = -2\\alpha_1\\Gamma_1(\\omega)^2sin(2\\pi\\omega) \\neq 0\n", "```\n", "What's $\\Gamma_1?$\n", "- for $\\alpha_1>0$ we have a maximum (plateau) at $\\omega=0$, i.e. red noise processes!\n", "```{figure} figures/L16/L16_5_spectrum_AR1.PNG\n", "---\n", "width: 30%\n", "---\n", "

Figure 5: Spectra of different AR(1)-processes.

\n", "```\n", "\n", "```{figure} figures/L16/L16_6_0.PNG\n", "---\n", "width: 30%\n", "---\n", "```\n", "```{figure} figures/L16/L16_6.PNG\n", "---\n", "width: 60%\n", "---\n", "

Figure 6: The spectra of AR(1)-processes with different $\\alpha_1$ compared with the spectrum of the driving white noise process. The spectra correspond to the time series atop.

\n", "```\n", "\n", "- fitting the AR(1)-process to a time series: AR(1)-processes or red noise are often chosen as a null hypothesis. AR(1)-processes are well defined by $\\sigma_{xy}$ and the lag-1 correlation because\n", "```{math}\n", "\\sigma_Z^2 = (1-\\alpha_1^2)\\sigma_{\\mathbf{X_t}}^2 = (1-\\rho_1^2)\\sigma_{\\mathbf{X_t}}^2\n", "```\n", "```{figure} figures/L16/L16_7_temperatures_series.PNG\n", "---\n", "width: 30%\n", "---\n", "```\n", "```{figure} figures/L16/L16_8_spectrum_temp.PNG\n", "---\n", "width: 60%\n", "---\n", "

Figure 7: Spectra of observed 24hrs mean temperature series (shown above). For comparison the fitted AR(1)-process spectra are shown together with the 95% confidence interval.

\n", "```\n", "\n", "```{figure} figures/L16/L16_9_ts_24h.PNG\n", "---\n", "width: 60%\n", "---\n", "

Figure 8: Left: Time series of 24hrs mean geopotential heights in Northern Germany (upper) and at the equatorial east Pacific (lower). Right: The auto-correlation functions corresponding to the time series in the left panels.

\n", "```\n", "```{figure} figures/L16/L16_10_spectrum_24h.PNG\n", "---\n", "width: 60%\n", "---\n", "

Figure 9: Spectra of observed 24h mean 500hPa geopotential heights time series shown in Fig. 8. For comparison the fitted AR(1)-process spectra are shown together with the 95% confidence interval.

\n", "```\n", "\n", "```{figure} figures/L16/L16_11_spectrum_AR1.PNG\n", "---\n", "width: 30%\n", "---\n", "

Figure 10: The spectrum of a discrete AR(1)-process with the time of one day averaged to a monthly mean time series (black, green). For comparison the spectrum resulting from eq.[9.5] of the AR(1)-rocess fitted to monthly mean time series is shown (red). What's eq.[9.5]?

\n", "```\n", "\n", "#### Spectrum of an AR(2)-process:\n", "- for $p=2$, [](#arp) can be rewritten as:\n", "```{math}\n", ":label: spec2\n", "\\Gamma(\\omega) = \\frac{\\sigma_Z^2}{1+\\alpha-1^2 + \\alpha_2^2 -2g(\\omega)},\\\\\n", "```\n", "- with:\n", "```{math}\n", "g(\\omega) = \\alpha_1(1-\\alpha_2)cos(2\\pi\\omega)+\\alpha_2cos(4\\pi\\omega)\n", "```\n", "```{figure} figures/L16/L16_12_0.PNG\n", "---\n", "width: 30%\n", "---\n", "```\n", "```{figure} figures/L16/L16_12_spectrum_AR2.PNG\n", "---\n", "width: 60%\n", "---\n", "

Figure 11: The spectra of AR(2)-processes with different $\\alpha_1$, $\\alpha_2$ (red line) and the estimated spectra from the time series above (black line), compared with the spectrum of the driving white process (green line) and a fitted AR(1)-process (dashed blue line).

\n", "```\n", "\n", "#### Estimating the spectra (periodogram)\n", "- let $\\{ \\mathbf{X_1}, \\mathbf{X_2}, ..., \\mathbf{X_T} \\}$ with t odd be a time series:\n", "```{math}\n", "\\mathbf{X_t} = A_0 + \\sum_{k=1}^{\\frac{T-1}{2}} \\left( a_kcos\\left(\\frac{2\\pi kt}{T}\\right) + b_ksin\\left(\\frac{2\\pi kt}{T}\\right) \\right)\n", "```\n", "- remember that th spectrum is a continous function of frequency and for the periodogram $a_k$ and $b_k$ are discrete\n", "```{math}\n", "\\begin{align*}\n", "A_0 &= \\hat{\\mu} = \\frac{1}{T} \\sum_{t=1}^T\\mathbf{X_t}\\\\\n", "a_j &= \\frac{2}{T} \\sum_{t=1}^T\\mathbf{X_t}cos(2\\pi \\omega_jt)\\\\\n", "b_j &= \\frac{2}{T} \\sum_{t=1}^T\\mathbf{X_t}sin(2\\pi \\omega_jt)\\\\\n", "\\omega_i &= \\frac{i}{T},~~~i=1,2,...,\\frac{T-1}{2}\n", "\\end{align*}\n", "```\n", "- the variance is given as:\n", "```{math}\n", "Var(\\mathbf{X_t}) = \\frac{1}{T} \\sum_{t=1}^T(\\mathbf{X_t}-\\mathbf{\\bar{X}})^2 = \\frac{1}{2} \\sum_{k=1}^{\\frac{T-1}{2}}(a_k^2+b_k^2)\n", "```\n", "- the elements $(a_k^2+b_k^2)$ are the periodogram of the finite time series $\\{ \\mathbf{X_1}, \\mathbf{X_2}, ..., \\mathbf{X_T} \\}$:\n", "```{math}\n", "I_{Tk} = \\frac{T}{4}(a_k^2+b_k^2)\n", "```\n", "- bad characteristics:\n", " 1. $I_{Tk} \\propto \\xi^2(2)$ relatively wide distribution, strongly skewed to larger values with peak at zero\n", " 2. Uncertainty of the spectrum coefficients is independent of the length of the time series\n", "```{figure} figures/L16/L16_13_spectrum_sine.PNG\n", "---\n", "width: 30%\n", "---\n", "

Figure 12: Spectral estimate of a sine function time series with 50 time steps. The red lines indicate uncertainties in the variance and period estimation.

\n", "```\n", "- best estimates:\n", " 1. Divide the time series into $m$ chunks of length $M=\\frac{T}{m}$.\n", " 2. Compute a periodogram $I_{Tk}^{(l)}$ for $k=0,1,..,\\frac{M}{2}$ from each chunk $l=1,2,..,m$\n", " 3. Estimate the spectrum by averaging the periodograms:\n", "```{math}\n", "\\hat{\\Gamma}(\\omega_j)=\\frac{1}{m}\\sum_{l=1}^m I_{Tk}^{(l)},\n", "```\n", "- estimator of the spectrum $\\propto \\xi^2(2m)$. estimate at each frequency is representative of a special bandwidth of $\\propto \\frac{1}{M}$\n", "\n", "```{figure} figures/L16/L16_14_periodogram.PNG\n", "---\n", "width: 60%\n", "---\n", "

Figure 13: Different estimates of the spectrum of AR(1)-process time series, based on the Periodogram.

\n", "```" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.6" } }, "nbformat": 4, "nbformat_minor": 5 }