<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/r-function-for-simulating-gaussian-processes/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://ntguardian.wordpress.com/2018/01/26/r-function-simulating-gaussian-processes/"> R – Curtis Miller's Personal Website</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>This semester my studies all involve one key mathematical object: <a rel="nofollow" href="https://en.wikipedia.org/wiki/Gaussian_process" target="_blank">Gaussian processes</a>. I’m taking a course on <a rel="nofollow" href="https://en.wikipedia.org/wiki/Stochastic_process" target="_blank">stochastic processes</a> (which will talk about <a rel="nofollow" href="https://en.wikipedia.org/wiki/Wiener_process" target="_blank">Wiener processes</a>, a type of Gaussian process and arguably the most common) and mathematical finance, which involves <a rel="nofollow" href="https://en.wikipedia.org/wiki/Stochastic_differential_equation" target="_blank">stochastic differential equations (SDEs)</a> used for derivative pricing, including in the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_model" target="_blank">Black-Scholes-Merton equation</a>. Then I’m involved in a Gaussian process and stochastic calculus reading group. So these processes will take up a lot of my attention.</p>
<p><span></span></p>
<p>In a conversation about these processes with a fellow graduate student I was explaining the idea that different kernels (covariance functions, or <img alt="E[X_t X_s]" title="E[X_t X_s]" src="https://s0.wp.com/latex.php?latex=E%5BX_t+X_s%5D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">) define different Gaussian processes and simply changing the kernel will produce new processes with completely different properties. Let <img alt="K(s, t)" title="K(s, t)" src="https://s0.wp.com/latex.php?latex=K%28s%2C+t%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> be the kernel of a process. <img alt="K(s, t) = min(s, t)" title="K(s, t) = min(s, t)" src="https://s0.wp.com/latex.php?latex=K%28s%2C+t%29+%3D+%5Cmin%28s%2C+t%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is the kernel associated with the Wiener process and produces a process that is continuous everywhere but not <a rel="nofollow" href="https://en.wikipedia.org/wiki/Differentiable_function" target="_blank">differentiable</a> anywhere, and with independent, Gaussian-distributed increments. On the other hand, the process defined by the kernel <img alt="K(s, t) = e^{-(s - t)^2}" title="K(s, t) = e^{-(s - t)^2}" src="https://s0.wp.com/latex.php?latex=K%28s%2C+t%29+%3D+e%5E%7B-%28s+-+t%29%5E2%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is not only continuous but differentiable everywhere, yet does not have independent increments.</p>
<p>I wanted to drive home the point that different kernels yield processes with wildly different properties by simulating and plotting them on a computer. So I whipped out the following R function in less than ten minutes (not counting documentation), and it does exactly what I want it to do.</p>
<pre>
library(MASS)

gaussprocess &lt;- function(from = 0, to = 1, K = function(s, t) {min(s, t)},
                         start = 0, m = 1000) {
  # Simulates a Gaussian process with a given kernel
  #
  # args:
  #   from: numeric for the starting location of the sequence
  #   to: numeric for the ending location of the sequence
  #   K: a function that corresponds to the kernel (covariance function) of
  #      the process; must give numeric outputs, and if this won't produce a
  #      positive semi-definite matrix, it could fail; default is a Wiener
  #      process
  #   start: numeric for the starting position of the process
  #   m: positive integer for the number of points in the process to simulate
  #
  # return:
  #   A data.frame with variables &quot;t&quot; for the time index and &quot;xt&quot; for the value
  #   of the process
  
  t &lt;- seq(from = from, to = to, length.out = m)
  Sigma &lt;- sapply(t, function(s1) {
    sapply(t, function(s2) {
      K(s1, s2)
    })
  })
  
  path &lt;- mvrnorm(mu = rep(0, times = m), Sigma = Sigma)
  path &lt;- path - path[1] + start  # Must always start at &quot;start&quot;
  
  return(data.frame(&quot;t&quot; = t, &quot;xt&quot; = path))
}
</pre>
<p>Below are example processes simulated by this function.</p>
<p><img alt="K(s, t) = min(s, t)" title="K(s, t) = min(s, t)" src="https://s0.wp.com/latex.php?latex=K%28s%2C+t%29+%3D+%5Cmin%28s%2C+t%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> (Wiener process)</p>
<p><img src="https://ntguardian.files.wordpress.com/2018/01/rplot07.png?w=456" data-recalc-dims="1"></p>
<p><img alt="K(s, t) = e^{-16(s - t)^2}" title="K(s, t) = e^{-16(s - t)^2}" src="https://s0.wp.com/latex.php?latex=K%28s%2C+t%29+%3D+e%5E%7B-16%28s+-+t%29%5E2%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> (Gaussian kernel)</p>
<p><img src="https://ntguardian.files.wordpress.com/2018/01/rplot08.png?w=456" data-recalc-dims="1"></p>
<p><img alt="K(s, t) = frac{1}{1 + left|s - t right|}" title="K(s, t) = frac{1}{1 + left|s - t right|}" src="https://s0.wp.com/latex.php?latex=K%28s%2C+t%29+%3D+%5Cfrac%7B1%7D%7B1+%2B+%5Cleft%7Cs+-+t+%5Cright%7C%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> (Something completely different)</p>
<p><img src="https://ntguardian.files.wordpress.com/2018/01/rplot09.png?w=456" data-recalc-dims="1"></p>
<p>Hopefully you found this code snippet entertaining, if not useful.</p>
<hr>
<p>I have created a video course published by <a rel="nofollow" href="http://packtpub.com/" target="_blank">Packt Publishing</a> entitled <a rel="nofollow" href="https://ntguardian.wordpress.com/video-courses/data-acquisition-and-manipulation-with-python/" target="_blank"><em>Data Acqusition and Manipulation with Python</em></a>, the second volume in a four-volume set of video courses entitled, <em>Taming Data with Python; Excelling as a Data Analyst</em>. This course covers more advanced <a rel="nofollow" href="http://pandas.pydata.org/" target="_blank">Pandas</a> topics such as reading in datasets in different formats and from databases, aggregation, and data wrangling. The course then transitions to cover getting data in “messy” formats from Web documents via web scraping. The course covers web scraping using <a rel="nofollow" href="https://www.crummy.com/software/BeautifulSoup/" target="_blank">BeautifulSoup</a>, <a rel="nofollow" href="http://www.seleniumhq.org/" target="_blank">Selenium</a>, and <a rel="nofollow" href="https://scrapy.org/" target="_blank">Scrapy</a>. If you are starting out using Python for data analysis or know someone who is, please consider <a rel="nofollow" href="https://www.packtpub.com/big-data-and-business-intelligence/data-acquisition-and-manipulation-python-video" target="_blank">buying my course</a> or at least spreading the word about it. You can buy the course directly or purchase a <a rel="nofollow" href="https://www.packtpub.com/mapt/" target="_blank">subscription to Mapt</a> and watch it there.</p>
<p>If you like my blog and would like to support it, spread the word (if not get a copy yourself)! Also, stay tuned for future courses I publish with Packt at the Video Courses section of my site.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://ntguardian.wordpress.com/2018/01/26/r-function-simulating-gaussian-processes/"> R – Curtis Miller's Personal Website</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/0Fzr02ip6qs" height="1"></a></div>