Module: Statsample::TimeSeries::Pacf

Included in:
Daru::Vector
Defined in:
lib/statsample-timeseries/timeseries/pacf.rb

Class Method Summary collapse

Class Method Details

.diag(mat) ⇒ Object

Returns diagonal elements of matrices



63
64
65
 
# 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) ⇒ Object

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



28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
 
# 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') ⇒ Object 



5
6
7
8
9
10
11
12
13
14
 
# 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) ⇒ Object

Solves matrix equations

Solves for X in AX = B



158
159
160
161
162
163
164
165
166
167
168
169
 
# 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) ⇒ Object

ToEplitz

Generates teoeplitz matrix from an array 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]]



135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
 
# 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') ⇒ Object

Yule Walker Algorithm

From the series, estimates AR(p)(autoregressive) parameter using Yule-Waler equation. See - 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



84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
 
# 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