function[mean_Ht,median_Ht,medianms_Ht,actual]=mean_median(data,nfor,S2,H,m,nh)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MATLAB code: 
% File: Mean_median.m
% Coded by: Dawit Zerom
%
% Returns the conditional mean and the conditional 
% median forecasts (Gaussian kernel); single- and multi-stage
%
% INPUT: 
%   nfor = Number of forecasts to generate
%   S2   = Size of the rolling forecast window
%   H    = Maximum forecast horizon
%   m    = Markov coefficient of the time series
%   nh   = Number of grid points for the bandwidth vector  
%
% OUTPUT: 
%   mean_Ht     = Single-stage conditional mean forecasts
%   median_Ht   = Single-stage conditional median forecasts
%   medianms_Ht = Multi-stage conditional median forecasts     
%   actual      = True values
%
% NOTE: The outputs are matrices each of size (nfor * H)
%
% References:
% De Gooijer, J.G., Gannoun, A., and Zerom, D. (2001).
%    Multi-stage kernel-based conditional quantile prediction in time series.
%    Communications in Statistics: Theory and Methods}, 30(12), 2499-2515.
%    DOI: 10.1081/sta-100108445.
% De Gooijer, J.G., Gannoun, A., and Zerom, D. (2002).
%    Mean squared error properties of the kernel-based multi-stage median 
%    predictor  for time series.
%    Statistics & Probability Letters, 56(1), 51-56.
%    DOI: 10.1016/S0167-7152(01)00169-9.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N1 = length(data);
N0 = N1 - nfor - H + 1;  % The first forecast origin

actual      = []; 
mean_Ht     = [];
median_Ht   = [];
medianms_Ht = [];

for t = N0-H:N1-2*H    % Start big loop
   % Bandwidth vector   
   hpre = 1.06*std(data(1:t)) * t^(-1/(m+1));
   hh   = linspace((0.1*hpre),(3*hpre),nh);   
   err2 = [];
    
   for i = 1:nh      
     if t == N0-H;
       for w = 1:S2
           temp = data(1:t-S2+w);
           for j = 1:H;
             [X,x,Y] = singledata(temp,j,m);             
             meany_H = meanfor(X,x,Y,hh(i));  % H-step ahead mean forecast            
             Ypred = data(t-S2+w+j) ;        
             err(w,j) = (Ypred - meany_H)^2;  % H-step ahead MSFE 
           end;                               % end loop j = 1:H
           meany_H = [];
       end;                                   % end loop w = 1:S2
        
       % Calculate the Euclidean norm to make the bandwidth a scalar    
         err1 = cumsum(err,2);
         [~,s2]=size(err1);
         errr  = err1(:,s2);
         err0  = errr;
         err1  = sum(err0);           
     elseif t >= N0 - H + 1;       
        temp = data(1:t);
        for j = 1:H
           [X,x,Y]   = singledata(temp,j,m);
           meany_H   = meanfor(X,x,Y,hh(i));
           Ypred     = data(t+j);  % -S2+w
           errm(1,j) = (Ypred - meany_H)^2;   
           meany_H   = [];
        end;                  % end loop j = 1:H    
        
        errer = cumsum(errm,2);
        errer = errer(H);
        s3    = length(err0);
        err0  = [err0(2:s3);errer];
        err1  = sum(err0);
     end;                       
     err2 = [err2;err1];     % An (nh * H) matrix
   end;                      % end loop i = 1:nh   
          
      [~,inx] = min(err2);
      hopk = hh(inx);  % Row vector (1 * H)      
      % The mean forecast at the optimal bandwidth
         temp      = data(1:t+H);         
%         [~,nn2] = size(temp);
         hopmed    = linspace(0.5*hopk,5*hopk,nh);      

    for jj=1:H 
        [X,x,Y]     = singledata(temp,jj,m);
        meany_H(jj) = meanfor(X,x,Y,hopk);
    end;               

  for ii = 1:nh;
     h = hopmed(ii);      
     for j = 1:H;
       [X,x,Y] = singledata(temp,j,m);
       % This is the optimal conditional mean j-step ahead forecast
       actualy(ii,j) = data(t+H+j);
       
       [nr,nc]   = size(X);
       [n01,n02] = size(x);      
       if ((n01==1)&&(n02==1));
            XX = ((X-x)/h ).^2; 
       else
            XX = ((X - vec2mat(x,nr))/h).^2;   
       end;      

       XX1 = cumsum(XX,2);
       A   = exp(-XX1(:,nc)/2)/(2*pi)^(m/2);              

%      Start median prediction         
       medd1   = meany_H(j);
       [nn,mm] = size(Y);
       if (mm > 1);
         Y2 = Y - vec2mat(medd1,nn);
       else
         Y2 = Y - medd1;
       end;
       YY    = abs(Y2);
       r     = mean(YY);
       toler = 100*r;
   
       while toler > 0.01*r;
         YY   = (YY<= 1e-4)*(1e-4)+YY.*((YY>1e-4));   
         knew = (1./YY).*vec2mat(A,mm);
         if sum(knew) == 0;
           k1 = 1e-10;
         else
           k1 = sum(knew);
         end;
 
         med10 = sum(Y.* knew)/k1;
         toler = abs(med10-medd1);
         medd1 = med10;
  
         [nn,mm] = size(Y);
         if (mm > 1);
           Y2 = Y - vec2mat(medd1,nn);
         else
           Y2 = Y - medd1 ;
         end;    
         YY = abs(Y2);   
       end;    
 
       med1       = medd1; 
       options(1) = -1;
       med1 = fmins('medianms1',med1,options,[],Y,h,A);
       med22(ii,j) = med1;
     end;  % end loop for  j = 1:H
  end;    % end loop for ii = 1:nh

   AA       = (med22 - actualy).^2;
   AA1      = cumsum(AA,2); AA1 = sqrt(AA1(:,H));
   [~,uu]   = min(AA1);
   alpha0   = med22(uu,:);
   actualy1 = actualy(1,:);
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
% Multihorizon prediction
   opmedmh     = median_mh(temp,meany_H,hopmed,H,m,actualy);

   actual      = [actual;actualy1];     
   mean_Ht     = [mean_Ht;meany_H];
   median_Ht   = [median_Ht;alpha0];
   medianms_Ht = [medianms_Ht;opmedmh];
end;  % end big loop for t = N0-H:N1-2*H      

%****************************************************************
% THE SUBFUNCTIONS        
function mean = meanfor(X,x,Y,h)
   [n1,n2]   = size(X);
   [n01,n02] = size(x);   
   if ((n01==1)&&(n02==1));
      XX = ((X-x)/h ).^2; 
   else
      XX = ((X - vec2mat(x,n1))/h).^2;   
   end;
   XX1 = cumsum(XX,2);
   A = exp(-XX1(:,n2)/2)/(2*pi)^(n2/2);     
   if sum(A)==0 ;
      A2 = 1e-300;
   else
      A2 = sum(A);
   end;  
   mean = sum(Y.* A)/A2;

function [Med,band] = median_mh(Z,p0,h0,H,m,act) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT: 
%   Z  = Time series up to the forecast origin of interest 
%   p0 = Optimal single-stage conditional Median forecast vector
%   h0 = Bandwidth vector from the single-stage forecast
%   H  = Maximum forecast horizon
%   m  = Markov coefficient of Z
%
% OUTPUT:  
%  med  = Optimal multi-stage Median prediction
%  band = Optimal bandwidth
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   n1 = length(Z);
   for u1 = 1:H;
      if (rem(n1-H+u1,H)==0);
         break;
      end;
   end;   
   u = u1;
   if (rem(m,H) == 0);
      q = m/H;
   else
      q = floor(m/H) + 1;
   end;    
   n = floor(n1/H);   
   for i = 1:n;
      c = i*H-u+1:1:i*H-u+H;
      W(i,:) = Z(c)';
   end;
   for i = 1:n-q+1;
     d       = i:1:i+q-1;
     xx      = (W(d,:))';
     X1(i,:) = xx(:)';
   end;   
   
   Y = W(1+q:1:n,:);
   X = X1(1:1:n-q,:);
   x = X1(n-q+1,:); 
   for j = 1:20;
      h1 = h0(j);   
      [n5,n6] = size(X);
      XX  = (X - vec2mat(x,n5))/h1;   
      XX1 = cumsum(XX.^2,2);
      A   = exp(-XX1(:,n6)/2)/(2*pi)^(q/2); 
      KK  = A;  
      medd2 = p0;  % This is the starting value

      [nn,mm] = size(Y);
      if (mm > 1);
        Y2 = Y - vec2mat(medd2,nn);
      else
        Y2 = Y - medd2;
      end; 
     
      YY = abs(Y2).^2;
      YY = sqrt(cumsum(YY,2));
      YY = YY(:,mm);
      r = norm(Y2);
      toler = r;
      ii = 1;

      while ((toler > (1e-3)*r) && (ii<=100));
        YY  = (YY <= 1e-4)*(1e-4) + YY.* ((YY>1e-4)|(YY<-1e-4));
        GRA = sum(vec2mat((1./YY),mm).*Y2.*vec2mat(KK,mm)); % Gradient
        ph1 = 0;
        for i = 1:nn;
          ph  = (1/YY(i))*(eye(mm)-(1/YY(i)^2)*Y2(i,:)*Y2(i,:)')*KK(i);  
          ph1 = ph1 + ph;
        end;
        
        HES     = ph1;     % Hessian matrix
        med20   = medd2 + GRA*pinv(HES);      
        toler   = norm(med20-medd2);
        medd2   = med20;
        [nn,mm] = size(Y);
        if (mm > 1);
          Y2 = Y - vec2mat(medd2,nn);
        else
          Y2 = Y - medd2;
        end;

        YY = sqrt(cumsum(Y2,2));
        YY = YY(:,mm);
        ii = ii+1;
      end;     
      med2 = medd2;
      options(1) = -1;
      med2 = fmins('medianms1',med2,options,[],Y,h1,KK);
      med221(j,:) = med2;
    end;   % end loop j = 1:20 

    AA  = (med221 - act).^2;
    AA1 = cumsum(AA,2); 
    AA1 = sqrt(AA1(:,H));
    [~,uu] = min(AA1);   
    Med      = med221(uu,:);
    band     = h0(uu);
    
function [X,x,Y] = singledata(Z,k,m)  
% INPUT: 
%   Z = Time series up to the forecast origin
%   k = a particular  forecast horizon
%   m = Markov coefficient
%
% OUTPUT:
%   X = Covariate matrix
%   Y = Response vector
%   x = Forecast origin covariate vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N  = length(Z);
n  = N-k-m+1;
X1 = Z(1:1:n);
 
for i=1:m-1;
    X1 = [X1 Z(1+i:1:n+i)];
end; 

X  = X1 ;
Y  = Z((1+m-1:1:n+m-1)+k);
x  = [X(n,2:m) Y(n)];        
    

function mat = vec2mat(vec,num) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VEC2MAT, Turns vectors into matrices (which are duplicates)
% INPUT:
% vec  = Row or column vector
% num  = Number of duplicates of vec to make
% OUTPUT:
% mat  = Matrix of duplicates of vec,
%          when vec is an nc * 1 column vector, mat is nc x num
%          when vec is a  1  * nr column vector, mat is num x nr
% NOTE:
%     To do this for BOTH a row and column vector, it is more 
%     convenient to use "meshgrid"
%     Copyright (c) J. S. Marron 1996, 1997
[nr,nc] = size(vec);
if nc == 1 && nr == 1;     %  Start with scalar, give warning
                           %  and return empty matrix
disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!');
disp(' Caution from vec2mat: not a vector, returning empty matrix');
disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!');
mat = [] ;
elseif nc == 1;     %  then treat as column vector
  mat = vec*ones(1,num);
elseif nr == 1;     %  then treat as row vector
  mat = ones(num,1)*vec;
else                %  then not a vector, give warning
                    %  and return empty matrix
disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!');
disp(' Caution from vec2mat: not a vector, returning empty matrix');
disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!');
  mat = [];
end; 

function [x, options] = fmins(funfcn,x,options,grad,varargin)
% FMINS Minimize function of several variables.
% FMINS has been replaced with FMINSEARCH. FMINS currently works but
% will be removed in the future. Use FMINSEARCH instead.
%
% X = FMINS('F',X0) attempts to return a vector X which is a local
% minimizer of F(x) near the starting vector X0. 'F' is a string
% containing the name of the objective function to be minimized,
% a string expression representing the objective function, or an
% inline function object. F(x) should be a scalar valued function
% of a vector variable.
%
% X = FMINS('F',X0,OPTIONS) uses a vector of control parameters. If
% OPTIONS(1) is positive, intermediate steps in the solution are
% displayed; the default is OPTIONS(1) = 0. OPTIONS(2) is the
% termination tolerance for x; the default is 1.e-4. OPTIONS(3) is
% the termination tolerance for F(x); the default is 1.e-4.
% OPTIONS(14) is the maximum number of function evaluations; the
% default is OPTIONS(14) = 200*length(x). The other components of
% OPTIONS are not used by FMINS. For more information, see FOPTIONS.
%
% X = FMINS('F',X0,OPTIONS,[],P1,P2,...) provides for additional
% arguments which are passed to the objective function, F(X,P1,P2,...)
% Pass an empty matrix for OPTIONS to use the default values.
%
% [X,OPTIONS] = FMINS(...) returns the number of function evaluations
% in OPTIONS(10).
%
% FMINS uses the Nelder-Mead simplex (direct search) method.
%
% See also FMIN, FOPTIONS.

% Reference: Jeffrey C. Lagarias, James A. Reeds, Margaret H. Wright,
% Paul E. Wright, "Convergence Properties of the Nelder-Mead Simplex
% Algorithm in Low Dimensions", May 1, 1997. To appear in the SIAM
% Journal of Optimization.
%
% Copyright 1984-2002 The MathWorks, Inc.
% $Revision: 5.23 $ $Date: 2002/04/08 20:26:45 $

if nargin<3, options = []; end
options = foptions(options);
prnt = options(1);
tolx = options(2);
tolf = options(3);
% The input argument grad is there for compatibility with FMINU in
% the Optimization Toolbox, but is not used by this function.

% Convert to inline function as needed.
funfcn = fcnchk(funfcn,length(varargin));

n = prod(size(x));
if (~options(14))
   options(14) = 200*n;
end;




