Using GP for autocorrelation in a time-series model: memory pressure, GP kernel definition

Modeling
, ,
14

functions {
matrix approx_L(int M, real scale, vector x, real l) {
int N = rows(x);

    real epsilon = 1 / (sqrt2() * l);
    real alpha = 1 / scale;
    real beta = (1.0 + (2.0 * epsilon / alpha)^2)^0.25;
    real delta =  alpha* sqrt((beta^2 - 1.0) / 2.0);
    vector[N] Ht[M];
    matrix[N, M] Htt;
    vector[N] xp = alpha * beta * x;
    real fsq = epsilon^2 / (alpha^2 + delta^2 + epsilon^2);
    real f = sqrt(fsq);
    Ht[1] = sqrt(sqrt(alpha^2 / (alpha^2 + delta^2 + epsilon^2))) * sqrt(beta) * exp(-delta^2 * x .* x);
    Ht[2] = f * sqrt(2.0) * xp .* Ht[1];
    for(n in 3:M) {
      Ht[n] = f * sqrt(2.0 / (n - 1.0)) * xp .* Ht[n - 1] - fsq * sqrt((n - 2.0) / (n - 1.0)) * Ht[n - 2];
    }

    for(n in 1:M) {
      Htt[:, n] = Ht[n];
    }
    return  Htt;
  }
}
data {
  int<lower=1> N; // number of (time) points in series
  vector[N] currents; // vector of Alongshore current speeds
  vector[N] wind; // vector of wind speeds
  //matrix[N,N] Dmat; // Distance matrix (scaled = /3600)

  int scale;  // set to 0.2
  vector[N] Distance;
  int<lower=1> M;  // approx points

  //prediction
  int<lower=1> S; // length of wind series to simulate over
  vector[S] sim_wind; // vector of wind values to simulate over
}
parameters {
  vector[M] g; // g is the variable intercept offset
  real a; // overall model intercept
  real b; // slope describing relationship to wind predictor
  real<lower=0> sigma_currents;
  real<lower=0> eta_sq;
  real<lower=0> inv_rho_sq;
  real<lower=0> sigma_sq;
}

transformed parameters {
  real<lower=0> rho_sq = inv(inv_rho_sq);
}

model {
    vector[N] mu;
    //priors
    eta_sq ~ cauchy(0,1);
    inv_rho_sq ~ cauchy(0,1);
    sigma_sq ~ cauchy(0,1);
    a ~ normal(0, 5);
    b ~ normal(0, 1);
    sigma_currents ~ normal(0, 10); // overall model sigma
    //model
    g ~ normal(0, 1);
    mu = a + approx_L(M, scale, Distance, inv_rho_sq) * g + b * wind;
    currents ~ normal(mu, sigma_currents);
}
generated quantities{
  vector[S] sim_cur; //simulated currents output
  sim_cur = a + b*sim_wind;
}

You can try this one. I modified the approx_L in, the approx_L is
from bbales2.

Choose M to something reasonable, since your frequencies are low, M doesn’t need to be very high.
Nyquist criterion.

Added, the mean of Distance should be zero. Not done in my code.

The approximation will work much better if the mean of the independent variable x is zero! The approximate kernels are not translation invariant or stationary or however you think about it! Here’s a picture of the first ten eigenfunctions (scaled by the square root of their eigenvalues). The approximation will fail if there is data out near the boundaries (< ~-1.2 or > ~1.2) cause all the eigenfunctions are near zero there:

1 Like