% InflationALS_PCE_11_23.m 
% Revised for new ALS with Global 5.0, 1/4/24
clear
kmin = 1   % fits AR(p), p = k - 1
kmax = 5
pmax = 4
realrate = 0     % compute real rates from interest rates if non-zero
                 % and 10-year interest rate provided
projections = 1  % = 1 for 10-yr proections, CIs
globflag = 1     % = 1 for Global significance stats

% If using CPI:
% P = xlsread('CPIAUCSL-2.xls');
% CPI-U, 1982-4 = 100, SA, 1/47 - 6/21, 894 obs.
% 1967-2006 is only to 1 decimal place, so inf has spurious patterns but I'll
%   use it anyway.  1947-66 is to 2 decimal places, probably uses 1967 = 100.
% t0str = '6/21'
% IndexStr = 'CPI-U'   % or 'PCE'
% t1 = 1947   % for CPI or 1959 for PCE
% BaseYr = '1982-4 = 100'  % for CPI

% If using PCE:
P = xlsread('PCEPI-30.xls');
% % PCE Chain-Type Price Index, 2017 = 100, SA, 1/59 - 11/23, 779 obs on P.
% xlsread finds a square array of numbers in the xls file.  Range doesn't
%   work in "basic" mode.
t0str = '11/23'
IndexStr = 'PCE'
t1 = 1959      % for PCE or 1947 for CPI
BaseYr = '2017 = 100'  % for PCE after 9/23

size(P)     % should be nP x 2, includes garbled date column 
P = P(:,2); % nPx1
nP = size(P,1)
PnP = P(nP)
tP = t1 + (0:nP-1)/12';  % nP x 1 gives dates of P readings

% Plot P level:
tmax = tP(nP) + 1
figure
hold on
plot([1945 2030], zeros(1, 2), 'k-', 'linewidth', 1.5)
plot(tP, P, 'linewidth', 1.5)
xlim([t1-1 tmax])
title([IndexStr ' SA'])
ylabel(BaseYr)
grid on

% Calculate and plot inflation:
inf = 1200.*log(P(2:nP)./P(1:nP-1)); %  ~ nP-1 x 1
tinf = tP(2:nP);
ninf = nP-1    % inf and tinf ~ ninf x 1, start in Feb 47 or Feb 59
figure
hold on
plot([1945 2030], [0 0], 'k-', 'linewidth', 1)
plot(tinf, inf, 'linewidth', 1.5)
xlim([t1-1 tmax])
title([IndexStr 'Inflation: SA, continuously compounded'])
ylabel('Annualized Percent Inflation')
grid on

disp(' ')
disp('Monthly')
disp('Using AR basis')
nn = ninf-pmax
tt = tinf((pmax+1):ninf);
y = inf((pmax+1):ninf);
XX = ones(nn,pmax+1);
if kmax > 1
  for k = 2:kmax
    p = k-1;
    XX(:,k) = inf(pmax-p+1  : ninf-p);
  end
end
%

%                                 ?

for k = kmin:kmax
  disp(' ')
  disp('Monthly')
  nn = nn
  k = k
  p = k-1
  X = XX(:,1:k);
%
% function [bf, Pf, bs, Ps, u, s2eps, rho, serho, NSR, NSR025, ...
%   NLR975, NLR,logl, lr, Global, pGlobal, C, nT, Tglob] ,,,
%   = ALS(x, y, tol, progflag, globflag)
% Global6.1 added 3/11/24.  
% 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," 1/26/22, 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
%
% 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)
% NSR = noise/signal s.d. ratio = 1/sqrt(rho)
% 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.     
% 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
% Global is kx1 vector of test statistics for beta(:, k) = 0, all t, if
%     globflag =1, or [ ] if globflag = 0.  Chi^2(n-k) under null.
% pGlobal is kx1 vextor of p-values for Global.
% C ~ nxnxkxk.  C(t1,t2,j1,j2) is Cov(bs(t1,j1),bs(t2,j2)), used to compute
% Global
% nT = size(Tglob), = DOF of Global stats, = round((n-k+1)/(2*NSR))
% Tglob ~ nTx1, indices for dates used to compute Global from C, bs 
%
[bf, Pf, bs, Ps, u, s2eps, rho, serho, NSR, NSR025, ...
    NSR975, NLR, logl, lr, Global, pGlobal, C, nT, Tglob] ...
      = ALS(X, y, .001, 0, globflag);
rho = rho
NSR = NSR
NSR025 = NSR025
NSR975 = NSR975
NLR = NLR
LRgain = 1/NLR
'LR for rho = 0:'
lr = lr
logl = logl
logl1500 = logl + 1500
[h,pjbstat, jbstat, critval5pct] = jbtest(u(k+1:nn));
jbstat = jbstat
pjbstat = pjbstat
chi2jbstat = chi2cdf(jbstat, 2, 'upper')
if globflag & k>1
  disp('Global test for beta(t,j) = 0')
  k = k
  p = p
  disp('j Global pGlobal')
  [(1:k)' Global pGlobal]
end
SSU = u(k+1:nn)'*u(k+1:nn)
s2eps = s2eps
sebf = nan(nn,k);
sebs = nan(nn,k);
for t = 1:nn
  sebf(t,:) = sqrt(diag(Pf(:,:,t)));
  sebs(t,:) = sqrt(diag(Ps(:,:,t)));
end
figure
hold on
title([IndexStr ' Filter coefficients,  AR(' num2str(k-1) ')'])
plot(tt, bf, 'linewidth', 1.5)
xlim([t1-1 tmax])
if k == 1
  ylim([-6 14])
end
grid on
legend ('Intercept', 'AR(1)', 'AR(2)', 'AR(3)', 'AR(4)', 'autoupdate', 'off')
plot([1947 2030], [0 0], 'k-', 'linewidth', 1)

if k > 2
  sumlags = sum(bf(:, 2:k)')';
  units = [0; ones(k-1,1)];
  sesumlags = nan(nn,1);
  for i = 1:nn
    covlags = Pf(:,:,i);
    sesumlags(i) = sqrt(units'*covlags*units); % or sum(sum(covlags))
  end 
  figure
  hold on
  title(['Sum of lag coefficients, ' IndexStr ', AR(' num2str(k-1) ')'])
  plot(tt, sumlags, 'b-', 'linewidth', 1.5)
  plot(tt, [sumlags+1.96*sesumlags sumlags-1.96*sesumlags], 'b-', 'linewidth', .5)
  legend ('Sum', '95% CI', 'autoupdate', 'off')
  plot([1945 2030], [0 0], 'k-', 'linewidth', 1)
  grid on
  xlim([t1-1 tmax])
end

if k > 1
  figure
  hold on
  title([IndexStr ' AR(' num2str(k-1) ') filter coefficient'])
  plot(tt, bf(:,k), 'b-', 'linewidth', 1.5)
  plot(tt, [bf(:,k)+1.96*sebf(:,k) bf(:,k)-1.96*sebf(:,k)], 'b-', 'linewidth', .5)
  legend (['AR(' num2str(k-1) ')'], '95% CI', 'autoupdate', 'off')
  plot([1945 2030], [0 0], 'k-', 'linewidth', 1)
  grid on
  xlim([t1-1 tmax])
  %
  figure
  hold on
  title([IndexStr ' AR(' num2str(k-1) ') Filter and Smoother Coefficients'])
  plot(tt, bf(:,k), 'b-', 'linewidth', 1.5)
  plot(tt, [bf(:,k)+1.96*sebf(:,k) ], 'b-', 'linewidth', .5)
  plot(tt, bs(:,k), 'r-', 'linewidth', 2.0)
  plot(tt, [bs(:,k)+1.96*sebs(:,k) ], 'r-', 'linewidth', 1.0)
  plot(tt, [ bf(:,k)-1.96*sebf(:,k)], 'b-', 'linewidth', .5)
  plot(tt, [bs(:,k)-1.96*sebs(:,k)], 'r-', 'linewidth', 1.0)
  legend ('Filter', '95% CI', 'Smoother', '95% CI', 'autoupdate', 'off')
  plot([1945 2030], [0 0], 'k-', 'linewidth', 1)
  grid on
  xlim([t1-1 tmax])
  %
  figure
  hold on
  title([IndexStr ' AR(' num2str(k-1) ') Smoother Z-stats'])
  plot(tt, bs(:,k)./sebs(:,k), 'r-', 'linewidth', 1.5)
  legend ('Smoother Z-stat', 'autoupdate', 'off')
  plot([1945 2030], [0 0], 'k-', 'linewidth', 1)
  grid on
  xlim([t1-1 tmax])
  if globflag & k>1
    Tglob = Tglob
    bsTglob = bs(Tglob,k)'
    sebsTglob = sebs(Tglob,k)'
    figure
    hold on
    title([IndexStr ' AR(' num2str(k-1) ') Smoother Z-stats'])
    plot(tt, bs(:,k)./sebs(:,k), 'r-', 'linewidth', 1.5)
    plot(tt(Tglob'), bs(Tglob',k)./sebs(Tglob',k), 'b.', 'markersize', 20) 
    legend ('Smoother Z-stat','Values used for Global test', 'autoupdate', 'off')
    plot([1945 2030], [0 0], 'k-', 'linewidth', 1)
    grid on
    xlim([t1-1 tmax])
  end
end
if k > 1
  if k == 2
    pibar = bf(:,1)./(1-bf(:,2));
    nomean = find(bf(:,2)>=1);
  else
    pibar = bf(:,1)./(1-sumlags);
    nomean = find(sumlags>=1);
  end
  pibar(nomean) = nan;
  nnomean = size(nomean,1);
  if nnomean > 0
    nomeanprime = nomean'
    ttnomean = tt(nomean)
    tt1900 = tt(nomean)-1900
  end
  figure
  hold on
  title(['Long Run Predicted ' IndexStr ' Inflation, AR(' num2str(k-1) ')'])
  plot(tt, pibar, 'linewidth', 1.5)
  if nnomean > 0
    plot(tt(nomean), zeros(nnomean,1), 'm*', 'markersize', 10)
    legend ('LR inf', 'Undefined', 'autoupdate', 'off')
  end
  plot([1945 2030], zeros(1, 2), 'k-', 'linewidth', 1)
  xlim([t1-1 tmax])
  grid on
end

if realrate
% Estimated 10-yr real rate
GS10 = xlsread('GS10.xlsx', 'B12:B830');
% 10-yr constant M Treasury yield, 4/53 - 6/21 from FRED
disp('size(GS10)')
size(GS10)
disp('size(pibar')
size(pibar)
disp('size(pibar(71:889)')
size(pibar(71:889))

figure
hold on
title('10-Year nominal Treasury yield, AL estimate of expected inflation')
plot(tt(71:889), [GS10 pibar(71:889)], 'linewidth', 2)
legend('10-Year Treasury yield', 'Expected inflation')
xlim([t1-1 tmax])
grid on
figure
hold on
title ('ALS estimate of expected 10-year real rate')
plot(tt(71:889), GS10 - pibar(71:889), 'linewidth', 2)
plot([1900 2100], [0 0], 'k-', 'linewidth', 1)
xlim([t1-1 tmax])
grid on

end % if realrate


piOneMo = nan(nn,1);
if k > 1
  for i = k:nn
    piOneMo(i) = bf(i,1) + y(i:-1:i-p+1)'*bf(i,2:k)';
  end  
  figure
  hold on
  title(['One month and long run predicted ' IndexStr ' inflation, AR(' num2str(k-1) ')'])
  plot(tt, piOneMo, 'r-', 'linewidth', 1.5)
  plot(tt, pibar, 'linewidth', 1.5)
  if nnomean > 0
    plot(tt(nomean), zeros(nnomean,1), 'm*', 'markersize', 10)
    legend ('1 mo inf', 'LR inf', 'LR undefined', 'autoupdate', 'off')
  else
    legend ('1 mo inf', 'LR inf', 'autoupdate', 'off')
  end
  plot([1945 2030], zeros(1, 2), 'k-', 'linewidth', 1)
  xlim([t1-1 tmax])
  grid on
  piOneMonn = piOneMo(nn)
  piLRnn = pibar(nn)
  
  eps = zeros(nn, 1);
  eps(k+1:nn) = inf(k+1:nn) - piOneMo(k:nn-1);
  cumunant = cumsum(eps)/12;    % Cumulative unant inf, %
  figure 
  hold on
  title (['Cumulative unanticipated ' IndexStr ' inflation, AR(' num2str(k-1) ')'])
  plot(tt, cumunant, 'b-', 'linewidth', 1.5)
  xlim([t1-1 tmax])
  grid on
  
  disp('One Month Predicted, LR Predicted, Cumulative Unanticipated, %')
  [tt' piOneMo pibar cumunant]'
  
end

  
if k == 1
  ttnn = tt(nn)
  bfttnn = bf(nn,1)
%   tt594 = tt(594)  % Nov. 08 cpi
%   bf594 = bf(594,1)
else
  ttnn = tt(nn)
  pibarttnn = pibar(nn)
%   tt594 = tt(594)
%   pibar594 = pibar(594)
end

figure
hold on 
title([IndexStr ' Scaled errors, AR(' num2str(k-1) ')'])
plot([1950 2030], zeros(1, 2), 'k-', 'linewidth', 1)
plot(tt, u, 'linewidth', 1.5)
xlim([t1-1 tmax])
grid on 

figure
hold on 
seps = sqrt(s2eps);
ustd = u/seps;
title([IndexStr ' Standardized scaled errors, AR(' num2str(k-1) ')'])
plot([1950 2030], zeros(1, 2), 'k-', 'linewidth', 1)
plot(tt, ustd, 'linewidth', 1.5)
xlim([t1-1 tmax])
grid on 


% % Projections.m
% function [yhat, ybar] = Projections(bf, pf, y0, rho, t0, t0str, M, reps, dt);
% % projects y(t+1) = b1 + b2 y(t) + ... for M periods of length dt
% % t0 = calendar date of last observation
% % t0str = string for t0
% % bf ~ 1xk
% % pf ~ k x k  = Cov(b)
% % y0 ~ k-1 x 1 = y(t), y(t-1), ... y(t-p) for p = k-1 
if k > 1 & projections > 0
  disp('Projections for ')
  k = k
  p = p
  i0 = nn
  [yhat, ybar, ybarMC1CIs, ybarMC2CIs] = ...
    Projections(bf(i0, :), Pf(:,:, i0), y(i0:-1:i0-p+1), ...
    rho, tt(i0), t0str, 120, 1000, 1/12);
end % if k > 1 & projections > 0
end % for k = 1:kmax

return





