|
{
A Handbook of Time Series Analysis, Signal Processing and Dynamics
D.S.G. POLLOCK
Queen Mary and Westfield College
The University of London
UK (pollok1.pdf)
- Delphi & Lazarus compatible (I changed the code to implement it correctly) -
}
type
{ TForm1 }
TA = Array of Array of Double;
procedure TForm1.LUSolve(start, n : integer; a : TA; x, b : Array of Double);
var
v, w, pivot, lambda : real;
i, j, k, pivotRow, g, h, finish : integer;
p : Array of integer;
d : Array of Double;
begin {LUSolveg}
finish := start + n - 1;
for i := start to finish do
begin {i; determine the scale factors}
p[i] := i;
d[i] := 0;
for j := start to finish do
if d[i] < Abs(a[i, j]) then
d[i] := Abs(a[i,j]);
end;
for i := start to finish - 1 do
begin {i; begin the process of reduction}
pivot := a[p[i], i];
for k := i + 1 to finish do
begin {k; search for a better pivot}
v := Abs(pivot)/d[p[i]];
w := Abs(a[p[k], i])/d[p[k]];
if v < w then
begin {interchange rows if better pivot is found}
pivot := a[p[k], i];
pivotRow := p[k];
p[k] := p[i];
p[i] := pivotRow;
end; {end interchange}
end; {k; end the search for a pivot}
for k := i + 1 to finish do
begin {k; eliminate a[k, i]}
lambda := a[p[k], i]/pivot;
for j := i + 1 to finish do
a[p[k], j] := a[p[k], j] - lambda * a[p[i], j];
a[p[k], i] := lambda; {save the multiplier}
end; {k}
end; {i; reduction completed}
for i := start to finish do
begin {i; forward-substitution}
x[i] := b[p[i]];
for j := i - 1 downto start do
x[i] := x[i] - a[p[i], j] * x[j];
end; {i; forward-substitution}
for i := finish downto start do
begin {i; back-substitution}
for j := i + 1 to finish do
x[i] := x[i] - a[p[i], j] * x[j];
x[i] := x[i]/a[p[i], i];
end; {i; back-substitution}
end; {LUSolve}
procedure TForm1.ARMACovariances(alpha, mu : Array of Double; gamma : Array of Double; varEpsilon : real; lags, p, q : integer);
var
i, j, t, tau, r : integer;
delta, f : Array of Double;
a : Array of Array of Double;
{Find the delta values}
begin
r := Max(p, q);
for t := 0 to q do
begin
delta[t] := mu[t] * varEpsilon;
tau := Min(t, p);
for i := 1 to tau do
delta[t] := delta[t] - delta[t - i] * alpha[i];
delta[t] := delta[t]/alpha[0]
end;
for i := 0 to r do
for j := 0 to r do
begin {i, j : form the matrix of alpha values}
a[i, j] := 0;
if ((i - j) >= 0) and ((i - j) <= p) then
a[i, j] := alpha[i - j];
if ((i + j) <= p) and (j > 0) then
a[i, j] := a[i, j] + alpha[i + j];
end; {i, j}
for i := 0 to r do
begin {i : form the RHS vector}
f [i] := 0.0;
for j := i to q do
f [i] := f [i] + mu[j] * delta[j - i]
end; {i}
{Solve for the initial autocovariances}
LUsolve(0, r + 1, a, gamma, f );
{Find the succeeding autocovariances}
for i := r + 1 to lags do
begin {i}
gamma[i] := 0.0;
for j := 1 to p do
gamma[i] := gamma[i] - alpha[j] * gamma[i - j];
gamma[i] := gamma[i]/alpha[0];
end; {i}
end; {ARMACovariances}
procedure TForm1.GaussNewtonARMA(p, q, n : integer; y : Array of Double; alpha, mu : Array of Double);
var
rhVec, delta, theta, newtheta, newAlpha, newMu : Array of Double;
crossYY, crossEY, crossEE : Array of Double;
storeY, storeE : Array of Double;
moments : TA;
stepScalar, oldSS, newSS : real;
i, iterations : integer;
convergence : boolean;
begin {GaussNewtonARMA}
convergence := false;
iterations := 0;
alpha[0] := 1;
mu[0] := 1;
newAlpha[0] := 1;
newMu[0] := 1;
{Prepare vectors using initial parameters}
for i := 0 to n do
storeE[i] := -y[i];
RationalExpansion(mu, q, n, n, storeE);
Convolution(alpha, storeE, p, n);
Covariances(storeE, storeE, crossEE, n, 0, 0);
oldSS := crossEE[0];
while (convergence = false) and (iterations < 20) do
begin {while : beginning of the major loop}
{Form the moment matrix and the RHS vector}
if q > 0 then
RationalExpansion(mu, q, n, n, storeE);
for i := 0 to n do
storeY[i] := y[i];
if q > 0 then
RationalExpansion(mu, q, n, n, storeY);
Covariances(storeY, storeY, crossYY, n, 0, p);
if q > 0 then
begin
Covariances(storeE, storeY, crossEY, n, p, q);
Covariances(storeE, storeE, crossEE, n, 0, q);
end;
MomentMatrix(crossYY, crossEE, crossEY, p, q, moments);
RHSVector(moments, crossYY, crossEY, alpha, p, q, rhVec);
{Find the updating vector}
Cholesky(p + q, moments, delta, rhVec);
{Update the value of theta}
stepScalar := -Sqrt(2);
repeat {until newSS < oldSS}
stepScalar := -stepScalar/sqrt(2);
for i := 1 to p + q do
delta[i] := stepScalar * delta[i];
for i := 1 to p do
newAlpha[i] := alpha[i] - delta[i];
for i := 1 to q do
newMu[i] := mu[i] - delta[p + i];
for i := 0 to n do
storeE[i] := -y[i];
if q > 0 then
RationalExpansion(newMu, q, n, n, storeE);
if p > 0 then
Convolution(newAlpha, storeE, p, n);
Covariances(storeE, storeE, crossEE, n, 0, 0);
newSS := crossEE[0];
until (newSS < oldSS * 1.0001);
iterations := iterations + 1;
oldSS := newSS;
for i := 1 to p + q do
begin {i}
if i <= p then
begin
alpha[i] := newAlpha[i];
theta[i] := alpha[i];
end
else if i > p then
begin
mu[i - p] := newMu[i - p];
theta[i] := mu[i - p];
end;
end; {i}
{Check for convergence}
if q = 0 then
convergence := true
else
convergence := CheckDelta(p + q, 0, delta, theta);
end; {while : end of the major loop}
end; {GaussNewtonARMA}
function CheckDelta(tolerance : real; q : integer; delta, mu : Array of Double) : boolean;
var
i, j : integer;
muNorm, deltaNorm : real;
begin
muNorm := 0;
deltaNorm := 0;
for i := 0 to q do
begin
muNorm := muNorm + Sqr(mu[i]);
deltaNorm := deltaNorm + Sqr(delta[i])
end;
if (deltaNorm/muNorm) > tolerance then
CheckDelta := false
else
CheckDelta := true;
end;{CheckDelta}
procedure TForm1.RationalExpansion(alpha : Array of Double; p, k, n : integer; beta : Array of Double);
var
i, j, r : integer;
begin
for j := 0 to n do
begin {j}
r := Min(p, j);
if j > k then beta[j] := 0;
for i := 1 to r do
beta[j] := beta[j] - alpha[i] * beta[j - i];
beta[j] := beta[j]/alpha[0]
end; {j}
end; {RationalExpansion}
procedure TForm1.Convolution(alpha, beta : Array of Double; p, k : integer);
var
gamma : real;
i, j, r, s : integer;
begin
for j := p + k downto 0 do
begin {j}
s := Min(j, p);
r := Max(0, j - k);
gamma := 0;
for i := r to s do
gamma := gamma + alpha[i] * beta[j - i];
beta[j] := gamma;
end; {j}
end; {Convolution}
procedure TForm1.Covariances(x, y : Array of Double; covar : Array of Double; n, p, q : integer);
var
t, j, s, f : integer;
begin {Covariances}
for j := -p to q do
begin {j}
s := Max(0, j);
f := Min(n, n + j);
covar[j] := 0;
for t := s to f do
covar[j] := covar[j] + x[t - j] * y[t];
covar[j] := covar[j]/(n + 1);
end; {j}
end; {Covariances}
procedure TForm1.MomentMatrix(covarYY, covarXX, covarXY : Array of Double; p, q : integer; moments : TA);
var
i, j : integer;
begin {MomentMatrix}
for i := 1 to p + q do
for j := 1 to p + q do
if i >= j then
begin {i; j : fill the lower triangle}
if (i <= p) and (j <= p) then
moments[i, j] := covarYY[i - j];
if (i > p) and (j <= p) then
moments[i, j] := covarXY[(i - p) - j];
if (i > p) and (j > p) then
moments[i, j] := covarXX[i - j];
moments[j, i] := moments[i, j]
end; {i, j}
end; {MomentMatrix}
procedure TForm1.RHSVector(moments : TA; covarYY, covarXY : Array of Double; alpha : Array of Double; p, q : integer; rhVec : Array of Double);
var
i, j : integer;
begin {RHSVector}
for i := 1 to p do
rhVec[i] := covarYY[i];
for i := p + 1 to p + q do
rhVec[i] := covarXY[i - p];
for i := 1 to p + q do
for j := 1 to p do
rhVec[i] := rhVec[i] + moments[i, j] * alpha[j];
end; {RHSVector}
procedure TForm1.Cholesky(n : integer; a : TA; x, b : Array of Double);
var
l : real;
i, j, k : integer;
begin {Cholesky}
for j := 1 to n do
for i := j to n do
begin {i; find the jth column of L}
l := a[i, j];
for k := 1 to j - 1 do
l := l - a[i, k] * a[j, k];
if i = j then
a[i, j] := Sqrt(l)
else
a[i, j] := l/a[j, j];
end; {i}
for i := 1 to n do
begin {i; forward-substitution}
x[i] := b[i];
for j := i - 1 downto 1 do
x[i] := x[i] - a[i, j] * x[j];
x[i] := x[i]/a[i, i];
end; {i; forward-substitution}
for i := n downto 1 do
begin {i; back-substitution}
for j := i + 1 to n do
x[i] := x[i] - a[j, i] * x[j];
x[i] := x[i]/a[i, i];
end; {i; back-substitution}
end; {Cholesky}
procedure TForm1.GramSchmidtPrediction(gamma : Array of Double; y : Array of Double; mu : TA; n, q : integer);
var
t : integer;
procedure MuLine(t : integer);
var
i, k : integer;
begin {t}
for i := 0 to t do
begin {i}
mu[t - i, t] := gamma[t - i];
for k := 0 to i - 1 do
mu[t - i, t] := mu[t - i, t]-mu[i - k, i] * mu[t - k, t] * mu[0, k];
if i < t then
mu[t - i, t] := mu[t - i, t]/mu[0, i];
end; {i}
end; {MuLine}
procedure Shiftmu;
var
t, j : integer;
begin
for t := 0 to q - 1 do
for j := 0 to t do
mu[j, t] := mu[j, t + 1];
end; {Shiftmu}
procedure PredictionError(q : integer);
var
k : integer;
begin
for k := 1 to q do
y[t] := y[t] - y[t - k] * mu[q - k, q];
end; {PredictionError}
begin {Main Procedure}
for t := 0 to q do
begin {t}
MuLine(t);
PredictionError(t);
end; {t}
for t := q + 1 to n do
begin {t}
Shiftmu;
MuLine(q);
PredictionError(q);
end; {Main Procedure}
end; {GSPrediction}
|
|
