% ALS.m, alsrho.m, hillright.m, alsloglik.m
function [bf, Pf, bs, Ps, u, s2eps, rho, serho, NSR, NSR025, ...
    NSR975, NLR,logl, lr, Global, pGlobal, C, nT, Tglob] ... 
    = ALS(x, y, tol, progflag, globflag)
% Global 6.2, 7/20/24
% Global 6.1, 3/11/24
% Global 6.0 added 3/1/24, 2/29/24 version archived as ALS0224Global50.m
% Global5.0 added 1/2/24.  1/1/24 version archived as ALS2023.m
% Converted to Matlab from GAUSS 1/9/12
% Gridsearch added 5/31/21
%
% Computes Adaptive Least Squares filter and smoother 
%     for model y = x*beta(t) + eps 
% for details, see J. Huston McCulloch, "Adaptive Least Squares:  
%     Recursive LS with Constant Noise-to-Signal Ratio," 2/xx/24, at
%     <http://www.asc.ohio-state/mcculloch.2/papers/ALS/>  
% Finds signal/noise variance ratio rho by ML  
%
% Inputs:
% x is nXk regressor matrix.  First column = ones(n,1) 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
% globflag = 1 to compute global significance test Global, pGlobal, C
%
% Outputs:
% 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 kXkxn super-matrix of coefficient variance matrices.  
%     Pf(:,:,t)) is covariance of b(t).  First k-1 pages are NaNs 
% 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 MLE estimated variance of eps 
% rho is NLE estimated signal/noise variance ratio. 
% serho = se(rho) based on curvature of logL at rho-hat
% NSR = noise/signal s.d. ratio = 1/sqrt(rho) = avg. lag from first
%     forecast
% NSR025 = lower bound of 95% CI based on likelihood ratio
% NSR975 = upper bound of 95% CI based on likelihood ratio
% NLR is asymptotic effective sample size.  Asymp. gain = 1/NLR.   
% seNLR = approx. se(NLR)  
% logl is maximized log likelihood 
% lr is 2*(logl(estimated rho) - logl(rho = 0)) 
%     This is informative, but distribution is not chi-square. 
%     5% critical value for constant parameters ~~ 2.3 in LLM
% Global is kx1 vector of test statistics for beta(:, k) = 0, all t, if
%     globflag =1, or [ ] if globflag = 0.  Chi^2(nT) under null.
% pGlobal is kx1 vextor of p-values for Global.
% C ~ n x n x k x k.  C(t1, t2,j1,j2) is Cov(bs(t1,j1), bs(t2,j2)), 
%     used to compute Global at nT values in Tglob.
% nT = size(Tglob), = DOF of Global stats.
% Tglob ~ nTx1 = t values used to compute Global from bs and C.
%
disp(' ')
disp('Function ALS')
format shortg     % was compact
debug = 0
k = size(x,2);
n = size(x,1);
if n ~= size(y,1)
  disp('Fatal error: row dimension of x and y not same in ALS')
  rowsx = size(x,1)
  rowsy = size(y,1)
  return
end
if n < k+1;
  disp('Fatal error: n < k+1 in ALS')
  n = n
  k = k
  return
end
% Preliminary gridsearch for best value of rho and plot:
NSRvec = n.^(1:-.01:0);   % = n ... 1, geometrically. 
rhovec = 1./(NSRvec.^2);     % 1/n^2 ... 1
loglvec = nan(1,101);
for i = 1:101
   [bf, Pf, bs, Ps, u, s2eps, loglvec(i)] = alsrho(x,y,rhovec(i),0,0,0);  
end
loglmax = max(loglvec);
imax = find(loglvec == loglmax);
if imax == 101
  disp('imax == 101 in ALS')    % debug to grid on
  figure
  plot(NSRvec, loglvec, 'r-')
  xlabel('NSRvec')
  ylabel('loglvec')
  grid on
  rhohi = rhovec(101);   % == 1
  rholo = rhovec(100);
elseif imax == 1
  disp('imax == 1 in ALS')       % debug to grid on
  figure
  plot(NSRvec, loglvec, 'r-')
  xlabel('NSRvec')
  ylabel('loglvec')
  grid on
  rholo = rhovec(1);      % == 1/sqrt(n)
  rhohi = rhovec(2);
else
  rholo = rhovec(imax-1);
  rhohi = rhovec(imax+1);
end
% Hill climb for best value of rho:
[rho, logl, d2logl] = hillright(@(rho) alsloglik(rho,x,y), rholo, rhohi, tol, progflag);
% compute logl0 for rho = 0:
[bf, Pf, bs, Ps, u, s2eps, logl0] = alsrho(x,y,0,0,0,0);   % that's a zero in logl0 
% smooflag = 1 in final pass:
[bf, Pf, bs, Ps, u, s2eps, logl, Global,pGlobal,C,nT,Tglob]...
      = alsrho(x,y,rho,1,globflag,debug);
lr = 2*(logl - logl0);
NLR = .5 + sqrt(.25 + 1/rho);
serho = sqrt(-1/d2logl);      
NSR = 1/sqrt(rho); 
figure
hold on
title ('LogL vs noise/signal ratio')
semilogx(NSRvec, loglvec, 'b-', 'linewidth', 1.5)
semilogx([1 n], [logl - 1.92 logl - 1.92], 'r-', 'linewidth', 1.5)
semilogx([NSR NSR], [logl logl-1.92], 'g-', 'linewidth', 1.5)
xlabel('Noise/Signal s.d. ratio, log scale')
ylabel('Log Likelihood')
legend('Log L', 'LR = 3.84', 'MLE')
grid on
%
if (loglvec(1) < logl-1.92 & loglvec(101) < logl-1.92)
  ilo = max(find(loglvec<logl-1.92 & NSRvec>NSR));  % since NSRvec is in decreasing order
  NSR975 = NSRvec(ilo) + (NSRvec(ilo+1) - NSRvec(ilo)) * ...
          (logl-1.92 - loglvec(ilo)) / (loglvec(ilo+1) - loglvec(ilo));
  ihi = min(find(loglvec<logl-1.92 & NSRvec<NSR));
  if debug 
    ihi = ihi
  end
  NSR025 = NSRvec(ihi-1) + (NSRvec(ihi) - NSRvec(ihi-1)) * ...
          (logl-1.92 - loglvec(ihi-1)) / (loglvec(ihi) - loglvec(ihi-1));
else
  NSR025 = [];
  NSR975 = [];
end
return

end       % function ALS


function [logl] = alsloglik (rho, x, y)
% Computes logl from rho for proc ALS to use with hillright.    
[bf, Pf, bs, Ps, u, s2eps, logl] = alsrho(x,y, rho, 0, 0, 0); 
return
end       % function alsloglik


function [bf, Pf, bs, Ps, u, s2eps, logl, Global, pGlobal, C, nT, Tglob] = ... 
    alsrho(x, y, rho, smooflag, globflag, debug)
% Global 6.2 7/20/24.
% Global 6.1 added 3/11/24.
%
% Computes Adaptive LS filter and/or smoother given rho, 
%     for model y = x*beta(t) + eps 
% Inputs: 
% 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
% globflag = 1 to compute Global, pGlobal, = 0 else.
% Outputs: 
% 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 kXkxn super-matrix of coefficient variance matices.  
%     Pf(t,:,:) is covariance of b(t), for 
%     t = k, ... n.  First k-1 pages are "NaN" 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. Rest should be NID(0, s2eps) 
%    if model correct.   
% s2eps is estimated variance of eps              
% logl is log likelihood
% Global is kx1 vector of test stats for beta(t,k) = 0, all t >= k
% pGlobal is kx1 vector of p-values for Global.  (~ Chi^2(n-k)).
% C ~ n x n x k x k.  C(t1,t2,j1,j2) = cov(bs(t1,j1), bs(t2,j2)).
% nT ~ 1x1 is # of t values in Tglob, = DOF of Global stats.
% Tglob ~ nT x 1 is set of t values from which Global computed using bs and
%     C
n = size(x,1);
k = size(x,2);
if n ~= size(y,1)
  disp('rows(x) =/ rows(y)')
  disp('rows(x)  rows(y)')
  [n size(y,1)]
  return
end
if n < k+1
  disp('rows(x) < cols(x) + 1')
  disp('rows(x)   cols(x)')
  [n k]
  return 
end
bs = [];
Ps = [];
Global = [];
pGlobal = [];
C = [];
nT = [];
Tglob = [];
% Compute filter, bf, Pf:
Nt = nan(n,1);
z = nan(k,n);
W = nan(k,k,n);
Wi = nan(k,k,n);
bf = nan(n,k); 
Pf = nan(k,k,n);
u = nan(n,1);
Xy = nan(n,k);    
XX = nan(k,k,n);
for t=1:n
 Xy(t,:) = x(t,:) * y(t);    % ~ nxk
 XX(:,:,t) = x(t,:)'*x(t,:); % ~ kxkxn
end
% Initialize filter:
logl = 0;
Nt(1) = 1;
z(:,1) = Xy(1,:)';
W(:,:,1) = XX(:,:,1); 
if k == 1 
  Wi(:,:,1) = inv(W(:,:,1));
  bf(1,1) = z(1,1)*Wi(1,1,1);
end
for t = 2:n
  fade = 1 + rho*Nt(t-1); 
  Nt(t) = Nt(t-1)/fade + 1;
  z(:,t) = z(:,t-1)/fade + Xy(t,:)';
  W(:,:,t) = W(:,:,t-1)/fade + XX(:,:,t);
  if t >= k
    Wi(:,:,t) = inv(W(:,:,t));
    bf(t,:) = (Wi(:,:,t)*z(:,t))';
  end
  if t > k
    s = sqrt(fade*x(t,:)*Wi(:,:,t-1)*x(t,:)' + 1);
    logl = logl - log(s);
    u(t) = (y(t) - x(t,:)*bf(t-1,:)') / s;
  end
end
  s2eps = (u((k+1):n)'*u((k+1):n)) / (n-k);
  logl = logl - ((n-k)/2) * (1 + log(2*pi*s2eps));
Pf(:,:,:) = s2eps * Wi(:,:,:);
%
% Compute Smoother, bs, Ps:
if smooflag | globflag;
  bs = nan(n,k);
  Ps = nan(k,k,n);
  zstar = nan(k,n);
  Wstar = nan(k,k,n);
  for t = 1:k-1
    [U S VV] = svd(W(:,:,t));
    for j = 1:t
      S(j,j) = 1/S(j,j);
    end
    for j = t+1:k
      % bugcheck:
      if abs(S(j,j)) > 1.e-10
        disp('S(j,j) ~= 0 for t < j <= k in SVD')
        t = t
        j = j
        Wt = W(:,:,t)
        U=U
        S=S
      end
      S(j,j) = 0;
    end
    Wi(:,:,t) = U*S*U';
  end
  V = nan(k,k,n);
  for t = 2:n
    V(:,:,t) = rho*Nt(t-1) * Wi(:,:,t-1);
  end
  Ws(:,:,n) = W(:,:,n);
  bs(n,:) = bf(n,:);
  Ps(:,:,n) = Pf(:,:,n);
  zstar(:,n) = Xy(n,:)';
  Wstar(:,:,n) = XX(:,:,n);
  for t = n-1:-1:k
    IWVi = inv(eye(k) + Wstar(:,:,t+1) * V(:,:,t+1));
    zstar(:,t) = IWVi * zstar(:,t+1) + Xy(t,:)';
    Wstar(:,:,t) = IWVi * Wstar(:,:,t+1) + XX(:,:,t);
    zs = z(:,t) + IWVi * zstar(:,t+1);
    Ws = W(:,:,t) + IWVi * Wstar(:,:,t+1); 
    Wsi = inv(Ws);
    bs(t,:) = (Wsi * zs)';
    Ps(:,:,t) = s2eps * Wsi;
  end   % for t 
end     % if smooflag or globflag
%
% Global 6.2:
if globflag
  disp('Global 6.2')
  n = n
  k = k
  if k == 1
    disp('Global not computed for k == 1')
    return
  end
  Phi = [y; zeros(k-1+k*(n-k),1)];
  Psi = [zeros(k-1, k*(n-k+1)) eye(k-1)];
  XXX = zeros(n-k+1, k*(n-k+1));
  for t = 1:n-k+1
    XXX(t, k*(t-1)+1:k*t) = x(t+k-1,:);
  end
  Psi = [Psi; [XXX zeros((n-k+1), k-1)]]; 
  Psi = [Psi; [x(1:k-1, :) zeros(k-1, k*(n-k)) -eye(k-1)]];
  II = eye(k*(n-k+1));  
  II = -II(1:k*(n-k), :) + II(k+1:k*(n-k+1), :);
  Psi = [Psi; [II zeros(k*(n-k), k-1)]];
  Xi = nan(k-1, k-1);
  for t = 1:k-1
    for tt = t:k-1
      VV = zeros(k,k);
      for s = tt+1:k
        VV = VV + V(:,:,s);
        Xi(t,tt) = x(t,:) * VV * x(tt,:)';
        Xi(tt,t) = Xi(t,tt);
      end % for s
    end % for tt
  end % for t   
  Xi = Xi
  Omega = eye(n+k-1+k*(n-k));
  Omega(n+1:n+k-1, n+1:n+k-1) = Xi;
  for t = k+1:n
    Omega(n+k+k*(t-k-1): n+k+k*(t-k)-1, n+k+k*(t-k-1): n+k+k*(t-k)-1) = ...
      V(:,:,t);          % V ~ kxkxn, = nans for t = 1
  end
  POPI = inv(Psi'*(Omega\Psi)); 
  Theta = (POPI*Psi')*(Omega\Phi);
  bsG = reshape(Theta(1:k*(n-k+1)), k, n-k+1)'; % ~ n-k+1xk.  bs~ nxk
  bsG = [nan(k-1,k); bsG];                      % ~ nxk.
  disp('Does bsG == bs?');
  figure
  hold on
  plot((1:n)', bs, 'b--', 'linewidth', 2)
  plot((1:n)', bsG, 'r-')
  title('bs, b--; bsG r-')
  grid on
  disp('max bs - bsG')
  err = max(max(abs(bs(k:n,:)- bsG(k:n,:)))) 
  %%% = 1.1e-12 for k = 2, PCE
  %                                          
  % Assemble C:
  % drop zeta portion:
  C = s2eps*POPI(1:k*(n-k+1), 1:k*(n-k+1));
  % pad with nans for bs(1:k-1, :) portions:
  C = [nan(k*(k-1),k*n); [nan(k*(n-k+1),k*(k-1)) C]]; 
  C = reshape(C,k,n,k,n);      % ~ kxnxkxn
  C = permute(C, [2 4 1 3]);   % ~ nxnxkxk
  %         bs ~ nxk, C(t1, t2, j1, j2) = cov(bs(t1,j1), bs(t2, j2))
  NSR = sqrt(1/rho);
  disp('Does Vs == VG?')
  Vs = nan(n,k);
  VG = nan(n,k);
  for t = 1:n
    for j = 1:k
      Vs(t,j) = Ps(j,j,t);   % Ps ~ kxkxn
      VG(t,j) = C(t,t,j,j);
    end
  end
  figure
  hold on
  plot((1:n)', Vs, 'b--', 'linewidth', 2)
  plot((1:n)', VG, 'r-')
  title ('Vs b--, VG r-')
  grid on
  disp('max Vs - VG')
  err = max(max(abs(Vs(k:n,:) - VG(k:n,:))))  
  %%  = 2.6e-13 for k=2, PCE, but e-5 for k>2!                     
  % Compute Glogal, pGlobal with sampling interval ~~ 2*NSR and
  %    reduce C to nT x nT x k x k:
  nT = round((n-k+1)/(2*NSR))
  Tglob = k-1 + round((.5:1:nT)*(n-k+1)/nT)
  Global = nan(k,1);
  for j = 1:k
    Global(j) = bs(Tglob,j)' * (C(Tglob,Tglob,j,j)\bs(Tglob,j));  % bs ~ nxk
  end
  pGlobal = chi2cdf(Global, nT, 'upper');
  disp('j, Global, pGlobal:')
  [(1:k)' Global pGlobal]
end % if globflag
%
return
end      % function alsrho

function [xmax,ymax,d2y] = hillright(f, xlo, xhi, tol, progflag)
% maximizes y = f(x) by golden ratio search, x >= xlo 
% loosely after Numerical Recipes.  
% call with hillright (@(x) f(x, data)), xlo, ...) 
%  d2y is 2nd derivative at xmax
%      should be negative if xmax > xlo  
%      = NaN 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. 
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;
  disp('Maximization progress');
end
%
while ((y1<y4) & (y2<y4) & (y3<y4)) 
  x2 = x3; y2 = y3;
  x3 = x4; y3 = y4;
  x4 = x3 + (x2 - x1);
  y4 = f(x4);
  if progflag
    disp(['expand right to ' num2str(x4)]) 
  end
end
%
if progflag 
  disp('contracting')
end
%
while dy > tol
if (y1 >= y2) & (y1 >= y3) & (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) || ((y2 == y3) & (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);
end
dy = max([y1; y2; y3; y4]) - min([y1; y2; y3; y4]);
if progflag;
  disp('x1 x2 x3 x4')
  [x1 x2 x3 x4 ]
  disp('y1 y2 y3 y4 dy')
  [y1 y2 y3 y4 dy] 
end
end
%
if y1 > y2;
    d2y = nan;
    xmax = x1;
    ymax = y1;
    return   
elseif y2 > y3;  % y1 <= y2 > y3 > y4
    dx2 = x2-x1;
    dx3 = x3-x1;
    dy2 = y2-y1;
    dy3 = y3-y1;
    d2y = 2*(dx2*dy3-dx3*dy2)/(dx2*dx3*(dx3-dx2));
    xmax = x2;
    ymax = y2;
    return
else;                  % y3 >= y2 >= y1 
    dx2 = x3-x2;
    dx3 = x4-x2;
    dy2 = y3-y2;
    dy3 = y4-y2; 
    d2y = 2*(dx2*dy3-dx3*dy2)/(dx2*dx3*(dx3-dx2)); 
    xmax = x3;
    ymax = y3;
    return
end
return
end % function hillright



