function [L] = gen_seasonal_trigonometric_L(k,s,m)
% GEN_SEASONAL_TRIGONOMETRIC_L - Gen(erate) Seasonal Trigonometric 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.
%   m    : the number of trigonometric components to keep/consider.
% 
% Our assumed model is a locally constant polynomial model with trigonometric predictors given by 
% 
% z_{n+j} = \beta_0 + \sum_{i=1}^k beta_i \frac{j^i}{i!} + \sum_{i=1}^{m} \beta_{1i} \sin( f_i t) + \beta_{2i} \cos( f_i t ) + \varepsilon_{n+j} 
% 
% With f_i = 2 \pi i / s 
% 
% See page 153 of the book by Abraham
% 
% Examples: 
% 
% The call MATLAB call: 
% 
% L = gen_seasonal_trigonometric_L(0,12,1)
% 
% will generate the matrix corresponding to Example 1 in the middle of page 154.
%   
% The call MATLAB call: 
% 
% L = gen_seasonal_trigonometric_L(1,12,2)
% 
% will generate the matrix corresponding to Example 2 in the middle of page 154.
%   
% 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.
% 
%-----

% create the polynomial portion of L: 
% 
L_1 = zeros( k+1 );
for jj=1:(k+1), 
  for ii=jj:(k+1), 
    L_1(ii,jj) = 1/factorial( ii-jj ); 
  end
end

% create the trigonometric portion of L: 
%
f = 2*pi/s; 

L_2 = zeros( 2*m );
for ii=0:(m-1),
  iip1 = ii+1; 
  L_2i = [  
       cos( f * iip1 ), sin( f * iip1 ), 
      -sin( f * iip1 ), cos( f * iip1 ) 
         ];
  L_2(2*ii+1,2*ii+1) = L_2i(1,1);
  L_2(2*ii+2,2*ii+1) = L_2i(2,1);
  L_2(2*ii+1,2*ii+2) = L_2i(1,2);
  L_2(2*ii+2,2*ii+2) = L_2i(2,2);
end

L = [ 
    L_1                  , zeros( k+1, 2 * m );
    zeros( 2 * m, k + 1 ), L_2 
    ];