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;