module Statsample::TimeSeries::Pacf

Public Class Methods

diag(mat) click to toggle source

Returns diagonal elements of matrices

# File lib/statsample-timeseries/timeseries/pacf.rb, line 63
def diag(mat)
  return mat.each_with_index(:diagonal).map { |x, r, c| x }
end
levinson_durbin(series, nlags = 10, is_acovf = false) click to toggle source

Levinson-Durbin Algorithm

Parameters

  • series: timeseries, or a series of autocovariances

  • nlags: integer(default: 10): largest lag to include in recursion or order of the AR process

  • is_acovf: boolean(default: false): series is timeseries if it is false, else contains autocavariances

Returns:

  • sigma_v: estimate of the error variance

  • arcoefs: AR coefficients

  • pacf: pacf function

  • sigma: some function

# File lib/statsample-timeseries/timeseries/pacf.rb, line 28
def levinson_durbin(series, nlags = 10, is_acovf = false)
  if is_acovf
    series = series.map(&:to_f)
  else
    #nlags = order(k) of AR in this case
    series = series.acvf.map(&:to_f)[0..nlags]
  end
  #phi = Array.new((nlags+1), 0.0) { Array.new(nlags+1, 0.0) }
  order = nlags
  phi = Matrix.zero(nlags + 1)
  sig = Array.new(nlags+1)

  #setting initial point for recursion:
  phi[1,1] = series[1]/series[0]
  #phi[1][1] = series[1]/series[0]
  sig[1] = series[0] - phi[1, 1] * series[1]

  2.upto(order).each do |k|
    phi[k, k] = (series[k] - (Statsample::Vector.new(phi[1...k, k-1]) * series[1...k].reverse.to_ts).sum) / sig[k-1]
    #some serious refinement needed in above for matrix manipulation. Will do today
    1.upto(k-1).each do |j|
      phi[j, k] = phi[j, k-1] - phi[k, k] * phi[k-j, k-1]
    end
    sig[k] = sig[k-1] * (1-phi[k, k] ** 2)

  end
  sigma_v = sig[-1]
  arcoefs_delta = phi.column(phi.column_size - 1)
  arcoefs = arcoefs_delta[1..arcoefs_delta.size]
  pacf = diag(phi)
  pacf[0] = 1.0
  return [sigma_v, arcoefs, pacf, sig, phi]
end
pacf_yw(timeseries, max_lags, method = 'yw') click to toggle source
# File lib/statsample-timeseries/timeseries/pacf.rb, line 5
def pacf_yw(timeseries, max_lags, method = 'yw')
  #partial autocorrelation by yule walker equations.
  #Inspiration: StatsModels
  pacf = [1.0]
  arr = timeseries.to_a
  (1..max_lags).map do |i|
    pacf << yule_walker(arr, i, method)[0][-1]
  end
  pacf
end
solve_matrix(matrix, out_vector) click to toggle source

Solves matrix equations

Solves for X in AX = B
# File lib/statsample-timeseries/timeseries/pacf.rb, line 158
def solve_matrix(matrix, out_vector)
  solution_vector = Array.new(out_vector.size, 0)
  matrix = matrix.to_a
  k = 0
  matrix.each do |row|
    row.each_with_index do |element, i|
      solution_vector[k] += element * 1.0 * out_vector[i]
    end
    k += 1
  end
  solution_vector
end
toeplitz(arr) click to toggle source

ToEplitz

Generates teoeplitz matrix from an array 
http://en.wikipedia.org/wiki/Toeplitz_matrix.
Toeplitz matrix are equal when they are stored in row & column major

== Parameters
  • arr: array of integers;

== Usage

  arr = [0,1,2,3]
  Pacf.toeplitz(arr)

  #=>  [[0, 1, 2, 3],
  #=>   [1, 0, 1, 2],
  #=>   [2, 1, 0, 1],
  #=>   [3, 2, 1, 0]]
# File lib/statsample-timeseries/timeseries/pacf.rb, line 135
def toeplitz(arr)
  eplitz_matrix = Array.new(arr.size) { Array.new(arr.size) }

  0.upto(arr.size - 1) do |i|
    j = 0
    index = i
    while i >= 0 do
      eplitz_matrix[index][j] = arr[i]
      j += 1
      i -= 1
    end
    i = index + 1; k = 1
    while i < arr.size do
      eplitz_matrix[index][j] = arr[k]
      i += 1; j += 1; k += 1
    end
  end
  eplitz_matrix
end
yule_walker(ts, k = 1, method='yw') click to toggle source

Yule Walker Algorithm

From the series, estimates AR(p)(autoregressive) parameter using 
Yule-Waler equation. See -
http://en.wikipedia.org/wiki/Autoregressive_moving_average_model

== Parameters
  • ts: timeseries

  • k: order, default = 1

  • method: can be 'yw' or 'mle'. If 'yw' then it is unbiased, denominator is (n - k)

== Returns
  • rho: autoregressive coefficients

  • sigma: sigma parameter

# File lib/statsample-timeseries/timeseries/pacf.rb, line 84
def yule_walker(ts, k = 1, method='yw')
  n = ts.size
  mean = (ts.inject(:+) / n)
  ts = ts.map { |t| t - mean }

  if method == 'yw'
    #unbiased => denominator = (n - k)
    denom =->(k) { n - k }
  else
    #mle
    #denominator => (n)
    denom =->(k) { n }
  end
  r = Array.new(k + 1) { 0.0 }
  r[0] = ts.map { |x| x**2 }.inject(:+).to_f / denom.call(0).to_f

  1.upto(k) do |l|
    r[l] = (ts[0...-l].zip(ts[l...n])).map do |x|
      x.inject(:*)
    end.inject(:+).to_f / denom.call(l).to_f
  end

  r_R = toeplitz(r[0...-1])

  mat = Matrix.columns(r_R).inverse
  phi = solve_matrix(mat, r[1..r.size])
  phi_vector = phi
  r_vector = r[1..-1]
  sigma = r[0] - (r_vector.map.with_index {|e,i| e*phi_vector[i] }).inject(:+)
  return [phi, sigma]
end