function [L] = gen_seasonal_indicator_L(k,s)
% GEN_SEASONAL_INDICATOR_L - Gen(erate) Seasonal Indicator L matrix
%
% Input: 
%   k    : order of the local polynomial model (k=1 is a local linear polynomial; k=2 is a local quadratic polynomial)  
%   s    : total number of seasons or seasonal indicators.
% 
% Our assumed model is a locally constant polynomial model with seasonal indicators given by 
% 
% z_{n+j} = \beta_0 + \sum_{i=1}^k beta_i \frac{j^i}{i!} + \sum_{i=1}^{s-1} \delta_i \mathrm{IND}_{ji} + \varepsilon_{n+j} 
% 
% See page 148 of the book by Abraham
% 
% The call: 
% 
% L = gen_seasonal_indicator_L(4,2)
% 
% will generate the matrix at the top of page 149.
%   
% Written by:
% -- 
% John L. Weatherwax                2009-06-01
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

L_11 = zeros( k+1 );
for jj=1:(k+1), 
  for ii=jj:(k+1), 
    L_11(ii,jj) = 1/factorial( ii-jj ); 
  end
end

L_21      = zeros( s-1, k+1 ); 
L_21(1,1) = 1; 

L_22 = zeros( s - 1 ); 
for jj=1:(s-1), 
  L_22(1,jj) = -1; 
end
for ii=2:(s-1), 
  L_22(ii,ii-1) = 1;
end

L = [ 
    L_11, zeros( k+1, s-1 ) ; 
    L_21, L_22 
    ];