proc(12) = als(x, y, tol, progflag);
@ Computes Adaptive Least Squares filter and smoother
for model y = x*beta(t) + eps @
@ for details, see J. Huston McCulloch, "The Kalman Foundations of
Adaptive Least Squares," Sept. 2004, at
@
@ Finds signal/noise variance ratio rho by ML @
@ x is nXk regressor matrix. First column = ones for constant term @
@ y is nX1 dependent variable vector @
@ tol is tolerance for precision of log likelihood. Try .001 or less @
@ progflag = 1 to show ML progress onscreen, = 0 for nothing onscreen @
@ returns {bf, pf, bs, ps, u, s2eps, rho, serho, tinf, setinf, logl, lr} @
@ bf is nXk matrix whose t-th row is b(t)', the
filter estimate of beta(t) conditional on y(1) ... y(t),
for t = k, ... n. First k-1 rows are "missing" code. @
@ pf is nXk^2 super-matrix of coefficient variance matrices.
reshape(pf[t,.],k,k) is covariance of b(t), for
t = k, ... n. First k-1 rows are "missing" code. @
@ bs, ps are similarly formatted smoother estimate and covariances. @
@ u is nx1 vector of variance-equalized forecasting errors.
Should be NID(0,s2eps) if model correctly specified. @
@ s2eps is estimated variance of eps @
@ rho is estimated signal/noise variance ratio. @
@ serho = se(rho) @
@ tinf is asymptotic effective sample size. Asymp. gain = 1/tinf. @
@ setinf = se(tinf) @
@ logl is maximized log likelihood @
@ lr is 2*(logl(estimated rho) - logl(rho = 0))
This is informative, but distribution is not chi-square. @
local k, n, rhohi, bf, pf, bs, ps, s2eps, logl0, logl, lr, rho, u, tinf ;
local serho, setinf, d2logl;
declare matrix alxglob;
declare matrix alyglob;
k = cols(x);
n = rows(x);
if n /= rows(y);
format /rdn 10,0;
print "rows(x) =/ rows(y) in ALS";
print "rows(x) " rows(x) " rows(y) " rows(y);
stop;
endif;
if n < k+1;
format /rdn 10,0;
print "rows(x) < cols(x) + 1 in ALS";
print "rows(x) " rows(x) " cols(x) " cols(x);
stop;
endif;
rhohi = 1/rows(y); @ makes T about sqrt(n) @
alxglob = x;
alyglob = y;
{rho, logl, d2logl} = hillright(&alsloglik, 0, rhohi, tol, progflag);
{bf, pf, bs, ps, u, s2eps, logl0} = alsrho(x,y,0,0); @ that's a zero in logl0 @
{bf, pf, bs, ps, u, s2eps, logl} = alsrho(x,y,rho,1);
lr = 2*(logl - logl0);
tinf = .5 + sqrt(.25 + 1/rho);
serho = sqrt(-1/d2logl);
setinf = serho/((1+4/rho)^.5 * rho^2);
retp(bf, pf, bs, ps, u, s2eps, rho, serho, tinf, setinf, logl, lr);
endp;
proc(1) = alsloglik (rho);
@ Computes logl from rho for proc als to use with hillright.
Receives x and y from als as globals alxglob, alyglob. @
local bf, pf, bs, ps, s2eps, logl, u;
{bf, pf, bs, ps, u, s2eps, logl} = alsrho(alxglob, alyglob, rho, 0);
retp(logl);
endp;
proc(7) = alsrho(x, y, rho, smooflag);
@ Computes Adaptive Learning filter and/or smoother given rho,
for model y = x*beta(t) + eps @
@ x is nXk regressor matrix. First column = ones for constant term @
@ y is nX1 dependent variable vector @
@ rho is signal/noise variance ratio @
@ smooflag = 1 to compute smoother, = 0 else. @
@ returns {bf, pf, bs, ps, u, s2eps, logl} @
@ bf is nXk matrix whose t-th row is b(t), the
filter estimate of beta(t) conditional on y(1) ... y(t),
for t = k, ... n. First k-1 rows are "missing" code. @
@ pf is nXk^2 super-matrix of coefficient variance matices.
reshape(pf[t,.],k,k) is covariance of b(t), for
t = k, ... n. First k-1 rows are "missing" code. @
@ bs, ps are similarly formatted smoother estimate and covariances,
if smooflag = 1. Null vectors returned if smooflag = 0. @
@ u is nx1 vector of variance-equalized forecast errors.
First k values are "missing" code. Should be NID(0, s2eps)
if model correct. @
@ s2eps is estimated variance of eps @
@ logl is log likelihood @
local k, n, tt, z, w, logl, s2eps, bf, pf, bs, ps, fade, s, wi, missing, u;
local zt, wt, zstar, wstar, ws, wsi, wiwrtw, tta, t;
missing = {.};
k = cols(x);
n = rows(x);
if n /= rows(y);
format /rdn 10,0;
print "rows(x) =/ rows(y)";
print "rows(x) " rows(x) " rows(y) " rows(y);
stop;
endif;
if n < k+1;
format /rdn 10,0;
print "rows(x) < cols(x) + 1";
print "rows(x) " rows(x) " cols(x) " cols(x);
stop;
endif;
z = zeros(k,1);
zt = {};
w = zeros(k,k);
wt = {};
logl = 0;
s2eps = 0;
bf = missing * zeros(n,k);
pf = missing * zeros(n,k*k);
u = missing * zeros(n,1);
bs = {};
ps = {};
tt = 0;
tta = {};
for t(1,n,1);
fade = (1+rho*tt);
if t > k;
s = sqrt(fade*x[t,.]*wi*x[t,.]' + 1);
logl = logl - ln(s);
u[t] = (y[t] - x[t,.]*bf[t-1,.]')/s;
endif;
tt = tt/fade + 1;
tta = tta|tt;
z = z/fade + x[t,.]'*y[t];
zt = zt~z;
w = w/fade + x[t,.]'*x[t,.];
wt = wt~vecr(w);
if t >= k;
wi = inv(w);
bf[t,.] = (wi*z)';
pf[t,.]=vecr(wi)';
endif;
endfor;
s2eps = (u[k+1:n]' * u[k+1:n]) / (n-k);
logl = logl - ((n-k)/2) * (1 + ln(2*pi*s2eps));
pf = pf * s2eps;
if smooflag;
bs = missing*zeros(n,k);
ps = missing*zeros(n,k*k);
zstar = zeros(k,1);
wstar = zeros(k,k);
t = n;
do until t < k;
w = reshape(wt[.,t],k,k);
wiwrtw = w*inv(w + rho*tta[t]*wstar);
ws = w + wiwrtw*wstar;
wsi = inv(ws);
ps[t,.] = s2eps*vecr(wsi)';
bs[t,.] = (wsi*(zt[.,t] + wiwrtw*zstar))';
zstar = wiwrtw*zstar + x[t,.]'*y[t];
wstar = wiwrtw*wstar + x[t,.]'*x[t,.];
t = t - 1;
endo;
endif;
retp(bf, pf, bs, ps, u, s2eps, logl);
endp;
proc(3) = hillright(&f, xlo, xhi, tol, progflag);
@ maximizes f(x) by golden ratio search @
@ loosely after Numerical Recipes. @
@ returns xmax, f(xmax), d2y
d2y is 2nd derivative at xmax
should be negative if xmax > xlo
= missing code if xmax = xlo @
@ may expand to right above xhi but not to left below xlo @
@ runs until dy < tol (Different than hill, hillur!) @
@ progflag = 1 for onscreen progress, 0 else. @
local f:proc, a, x1, x2, x3, x4, dx, y1, y2, y3, y4, dy;
local dx2, dx3, dy2, dy3, d2y, missing;
a = (3 - sqrt(5))/2; @ a = .382... = 1 - golden ratio @
x1 = xlo;
x4 = xhi;
dx = x4 - x1;
x2 = x1 + a * dx;
x3 = x4 - a * dx;
y1 = f(x1); y2 = f(x2); y3 = f(x3); y4 = f(x4);
dy = tol + 1000;
if progflag;
format /rdn 13,6;
print;
"Maximization progress";
endif;
right:
if y1 < y4 and y2 < y4 and y3 < y4;
x2 = x3; y2 = y3;
x3 = x4; y3 = y4;
x4 = x3 + (x2 - x1);
y4 = f(x4);
if progflag; print "expand right to " x4 y4; endif;
goto right;
endif;
if progflag; print "contracting"; endif;
do until dy < tol;
if (y1 >= y2) and (y1 >= y3) and (y2 <= y3);
x4 = x2; y4 = y2;
x2 = x1 + a*(x4 - x1);
x3 = x4 - a*(x4 - x1);
y2 = f(x2); y3 = f(x3);
elseif (y2 > y3) or ((y2 == y3) and (y1 >= y4));
x4 = x3; y4 = y3;
x3 = x2; y3 = y2;
x2 = x1 + (x4 - x3);
y2 = f(x2);
else;
x1 = x2; y1 = y2;
x2 = x3; y2 = y3;
x3 = x4 - (x2 - x1);
y3 = f(x3);
endif;
dy = maxc(y1|y2|y3|y4) - minc(y1|y2|y3|y4);
if progflag;
print "x1 x2 x3 x4 dy" x1 x2 x3 x4 dy;
print "y1 y2 y3 y4 " y1 y2 y3 y4 ;
endif;
endo;
missing = {.};
if y1 > y2;
d2y = missing;
retp (x1, y1, d2y);
elseif y4 > y3; @ this case should never occur @
d2y = missing;
retp (x4, y4, d2y);
elseif y2 > y3;
dx2 = x2-x1;
dx3 = x3-x1;
dy2 = y2-y1;
dy3 = y3-y1;
d2y = 2*(dx2*dy3-dx3*dy2)/(dx2*dx3*(dx3-dx2));
retp (x2, y2, d2y);
else; @ y3 >= y2 @
dx2 = x3-x2;
dx3 = x4-x2;
dy2 = y3-y2;
dy3 = y4-y2;
d2y = 2*(dx2*dy3-dx3*dy2)/(dx2*dx3*(dx3-dx2));
retp (x3, y3, d2y);
endif;
endp;