% SiZer1d (Family and SiZer1d Plots for fixed gammahat(k) function and nstar approx.)
%
% GPSZ1D, General Purpose Significant derivative Zero crossings for dep. data
%     Creates gray level map (function of location and bandwidth),
%     which is:
%          1 (very dark)    - Deriv. Sig. > 0 
%          2 (darker gray)  - Eff. SS < 5
%          3 (lighter gray) - Eff. SS >= 5, but CI contains 0
%          4 (very light)   - Deriv. Sig. < 0 
%         
%

clear;

tic;


% data
%
load m5har2pt.txt;                                %%% to be considered input %%%
%load m5dt.txt;

ydat=m5har2pt;           % Yi's: column vector n*1
n=size(ydat,1);               % sample size n
samplesize=int2str(n);


model=m5har2pt;                                 % deterministic Yi's

xdat=(1:n)';                  % Xi's: column vector n*1
Z=[xdat ydat];                % matrix of data n*2


% bandwidth vector (very wide range of bandwidths)
%
hmin=2;
hmax=n-1;
nh=11;                                           %%% to be considered input %%%
vh=logspace(log10(hmin),log10(hmax),nh);


% to calculate vybar
%
nn=fix(sqrt(n));
p=n-nn^2;


% to compute the autocovariance function
%
kk=n-1;                % suggested maximum: lag fix(n/4), but it does not work 
maxlag=int2str(kk);


figure(1)

subplot(4,3,1)
gpnpr(Z,vh,[1;n;n]);
hold on;
plot(xdat,ydat,'.');
hold off;
title('Family Plot','FontSize',12);
axis([1 n min(ydat)-10 max(ydat)+10]);
set(gca,'XTick',[1 fix(n/2) n]);


% binned data
%
b=ones(size(xdat));
bincts=[b ydat];


% construction of surfaces
%
MHAT=[];          % matrix of estimated regression function (n*nh) 
MDHAT=[];         % matrix of estimated derivative regression function (n*nh)

ACOV=[];          % matrix of sample autocovariance function ((kk+1)*nh)
ACOR=[];          % matrix of sample autocorrelation function ((kk+1)*nh)         

KVEC0=[];         % n*nh
KVEC1=[];         % n*nh
DEN=[];           % n*nh
VNUM1=[];         % n*nh
S0=[];            % n*nh
S1=[];            % n*nh

ESS=[];           % n*nh

EHAT=[];          % n*nh
SUMEHAT2=[];      % 1*nh

VYBAR=[];         % 1*nh

NSTAR=[];         % 1*nh

% loop through bandwidths: 
% first loop to calculate what does not depend on gammahat
%
for ih=1:nh;
  h=vh(ih);

  % kernel vector
  %
  arg=linspace(0,(n-1)/h,n)';         % vector n*1 [(i/h), i=0,...,n-1] 
  argh=arg*h;                         % vector n*1 [i, i=0,...,n]  

  kvec0=exp(-(arg.^2)/2)/sqrt(2*pi);  % vector n*1 [K(i/h), i=0,...,n-1, K~N]
  kvec1=kvec0.*argh;                  % vector n*1 [K(i/h)*(i), i=0,...,n-1]
  kvec2=kvec1.*argh;                  % vector n*1 [K(i/h)*(i^2), i=0,...,n-1]
 
  skvec0=kvec0.^2;                    % vector n*1 [(K(i/h))^2, i=0,...,n-1]
  skvec1=skvec0.*argh;                % vector n*1 [(K(i/h))^2*(i), i=0,...,n-1]
  skvec2=skvec1.*argh;                % vector n*1 [(K(i/h))^2*(i^2),i=0,..,n-1]


  % s0,s1,s2,t0,t1,u0,u1,u2
  %
  kvec0=[flipud(kvec0(2:n));kvec0];     % vector (2*n-1)*1 of kernel values

  KVEC0=[KVEC0 kvec0];

  skvec0=[flipud(skvec0(2:n));skvec0];  % vector (2*n-1)*1 of square kernels

  s0=conv(bincts(:,1),kvec0);         % vector [3*n-2]*1 
  s0=s0(n:2*n-1);                     % vector n*1

  S0=[S0 s0]; 

  t0=conv(bincts(:,2),kvec0);         % vector [3*n-2]*1 
  t0=t0(n:2*n-1);                     % vector n*1 

  u0=conv(bincts(:,1),skvec0);        % vector [3*n-2]*1
  u0=u0(n:2*n-1);                     % vector n*1 

  kvec1=[-flipud(kvec1(2:n));kvec1];    % vector (2*n-1)*1 of kernel*argh

  KVEC1=[KVEC1 kvec1];

  skvec1=[-flipud(skvec1(2:n));skvec1]; % vector (2*n-1)*1 of kernel^2*argh

  s1=-conv(bincts(:,1),kvec1);        % vector [3*n-2]*1 
  s1=s1(n:2*n-1);                     % vector n*1 

  S1=[S1 s1]; 

  t1=-conv(bincts(:,2),kvec1);        % vector [3*n-2]*1 
  t1=t1(n:2*n-1);                     % vector n*1 

  u1=-conv(bincts(:,1),skvec1);       % vector [3*n-2]*1 
  u1=u1(n:2*n-1);                     % vector n*1

  kvec2=[flipud(kvec2(2:n));kvec2];     % vector (2*n-1)*1 of kernel*argh^2
  skvec2=[flipud(skvec2(2:n));skvec2];  % vector (2*n-1)*1 of kernel^2*argh^2

  s2=conv(bincts(:,1),kvec2);         % vector [3*n-2]*1
  s2=s2(n:2*n-1);                     % vector n*1 

  u2=conv(bincts(:,1),skvec2);        % vector [3*n-2]*1 
  u2=u2(n:2*n-1);                     % vector n*1


  % estimation of the regression function and its derivative
  %
  den=s2.*s0-s1.*s1; 

  DEN=[DEN den];

  mhat=(s2.*t0-s1.*t1)./den;          % vector n*1 of estimated values for m
  mdhat=(s0.*t1-s1.*t0)./den;         % vector n*1 of estimated values for m' 

  MHAT=[MHAT mhat];                   % it will become a matrix n*nh 
  MDHAT=[MDHAT mdhat];                % it will become a matrix n*nh 


  % residuals
  %
  ehat=ydat-mhat;
  EHAT=[EHAT ehat];

  
  % sum of the squared residuals
  %
  ehat2=ehat.^2;
  sumehat2=sum(ehat2);
  SUMEHAT2=[SUMEHAT2 sumehat2];

  
  % vybar
  %
  YBAR=[];
 
if p<nn
 
  for i=1:p;
    ybar=sum(ehat((nn+1)*(i-1)+1:(nn+1)*i))/(nn+1);
    YBAR=[YBAR ybar];
  end
  
  for i=p+1:nn;
    ybar=sum(ehat(nn*(i-1)+p+1:nn*i+p))/nn;
    YBAR=[YBAR ybar];
  end
  
  ybartot=sum(YBAR)/nn;
  
  vybar=sum((YBAR-ybartot).^2)*(1/(nn*(nn-1)));
  
  VYBAR=[VYBAR vybar];
  
else

  pn=p-nn;

  for i=1:pn;
    ybar=sum(ehat((nn+2)*(i-1)+1:(nn+2)*i))/(nn+2);
    YBAR=[YBAR ybar];
  end
  
  for i=pn+1:nn;
    ybar=sum(ehat((nn+1)*(i-1)+pn+1:(nn+1)*i+pn))/(nn+1);
    YBAR=[YBAR ybar];
  end
  
  ybartot=sum(YBAR)/nn;
  
  vybar=sum((YBAR-ybartot).^2)*(1/(nn*(nn-1)));
  
  VYBAR=[VYBAR vybar];

end
    
  % sample autocovariance and autocorrelation function for the residuals
  %  
  acov=autocov(ehat,kk);
  acov=acov';
  g0=acov(1);
  acor=acov/g0;

  ACOV=[ACOV acov];                  % it will become a matrix kk*nh
  ACOR=[ACOR acor];                  % it will become a matrix kk*nh
  
  
  % estimation of part of the numerator of the variance
  %
  vnum1=g0*((s0.^2).*u2+(s1.^2).*u0-2*s0.*s1.*u1);
  VNUM1=[VNUM1 vnum1];               % it will become a matrix n*nh


  % effective sample size without taking into account the dependence
  %
  ess=sqrt(2*pi)*s0;                 % n*1
  ESS=[ESS ess];                     % it will become a matrix n*nh


  % quantities for effective sample size taking into account the dependence
  %
  nstar=acov(1)/vybar;
  NSTAR=[NSTAR nstar];
  
end;


% loop through bandwidths: 
% second loop to calculate what does depend on gammahat and the corresponding
% SiZer1d Plot
%
for i=1:nh;

  VMDHAT=[];        % matrix of estimated variance of the derivative function 
  GQUANT=[];
  
  acov=ACOV(:,i);   % fixed gammahat function
  
  nstar=NSTAR(:,i);
  
  
  % effective sample size
  % 
  ESTAR=(nstar/n)*ESS; 
   
   
  % loop through bandwidths: 
  % continuation of the first loop, once gammahat is fixed
  %

  for j=1:nh;

    kvec0=KVEC0(:,j);
    kvec1=KVEC1(:,j);

    % estimation of the variance of the derivative of the regression function
    % 

    % w00, w01, w10, w11
    %
    ww00=[];
    ww11=[];
    ww01=[];
    ww10=[];

    SS0=[];
    SS1=[];

    kvec1f=flipud(kvec1);

    kkvec0=[0;kvec0(1:n-1)];                   % n*1
    kkvec1=[0;kvec1f(1:n-1)];                  % n*1

    KA=kkvec0*acov(3:kk+1)';                   % n*(kk-2)
    KA1=kkvec1*acov(3:kk+1)';                  % n*(kk-2)

    ss0=conv(acov(2:kk+1),kvec0(1:2*n-2));     % (3n-4)*1
    ss1=conv(acov(2:kk+1),kvec1f(1:2*n-2));    % (3n-4)*1


    for k=1:n-kk;

      ss00=ss0(n-k:2*n-(k+1));                
      ss11=ss1(n-k:2*n-(k+1));

      SS0=[SS0;ss00']; 
      SS1=[SS1;ss11'];

    end;   
 

    for k=n-kk+1:n-1; 

      ss00=[ss0(n-k);ss00(1:n-1)]-KA(:,n-k); 
      SS0=[SS0;ss00']; 

      ss11=[ss1(n-k);ss11(1:n-1)]-KA1(:,n-k); 
      SS1=[SS1;ss11']; 
    
    end;


    for k=1:n;

      kk0=kvec0(n-(k-1):2*n-(k+1));
      kk1=kvec1(n-(k-1):2*n-(k+1));

      w00=kk0'*SS0(:,k);
      w11=kk1'*SS1(:,k);
      w01=kk0'*SS1(:,k);
      w10=kk1'*SS0(:,k);

      ww00=[ww00;w00];
      ww11=[ww11;w11];
      ww01=[ww01;w01];
      ww10=[ww10;w10];

    end;


    den=DEN(:,j);
    vnum1=VNUM1(:,j);
    s0=S0(:,j);
    s1=S1(:,j);

    vnum2=2*((s0.^2).*ww11+(s1.^2).*ww00-s0.*s1.*(ww01+ww10));

    vmdhat=(vnum1+vnum2)./(den.^2); 

    VMDHAT=[VMDHAT vmdhat];       % it will become a matrix n*nh
 
 
    % quantiles for simultaneous CI
    %
    alpha=0.05;

    estar=ESTAR(:,j);

    flag=(estar>=5);
    if sum(flag)>0
     
      mestar=mean(estar(flag));
      m=n/mestar;                  % effective number of independent groups
       
      beta=(1-alpha)^(1/m);
    
      gquant=-phiinv((1-beta)/2);
      
    else
    
      gquant=inf;
      
    end;
    
    GQUANT=[GQUANT gquant];
    
  end;

  
  % construction of the scale space CI surface
  %

  LOWCI=MDHAT-vec2mat(GQUANT,n).*sqrt(VMDHAT);   % lower confidence surface
  UPCI=MDHAT+vec2mat(GQUANT,n).*sqrt(VMDHAT);    % upper confidence surface

  % construction "gray level map":
  %          1 (very dark)    - Deriv. Sig. > 0 
  %          2 (darker gray)  - Eff. SS < 5
  %          3 (lighter gray) - Eff. SS >= 5, but CI contains 0
  %          4 (very light)   - Deriv. Sig. < 0 

	  

  MAPOUT=[];

  MAPOUT=3*ones(size(LOWCI));            % default is lighter grey

  flag=(LOWCI>0);
  ssflag=sum(sum(flag));
  if ssflag>0;
    MAPOUT(flag)=ones(ssflag,1);         % dark grey where significantly >0
  end;

  flag=(UPCI<0);
  ssflag=sum(sum(flag));
  if ssflag>0;
    MAPOUT(flag)=4*ones(ssflag,1);       % light grey where significantly <0
  end;

  flag = (ESTAR<=5) ;
  ssflag = sum(sum(flag)) ;
  if ssflag>0 ;
    MAPOUT(flag) =2*ones(ssflag,1) ;     % middle darker grey where estar <5
  end ;

  
  MAPOUT=MAPOUT';
  MAPOUT=flipud(MAPOUT);

  comap=[0,0,1;.35, .35, .35;.5,0,.5;1,0,0];

  subplot(4,3,i+1)

  tt=int2str(i);

  %image([1,n],[log10(hmax),log10(hmin)],MAPOUT);
  image(MAPOUT);
  colormap(comap);
  ylabel('log10(h)');
  title(['hp(' tt ')'],'FontSize',12)
  set(gca,'XTick',[1 fix(n/2) n])
  set(gca,'YTick',[])
  

end;


  orient tall; 
%  print -dpsc c:\utenti\vita\articoli\cfig4
print -dpsc CFIG4.ps
print -depsc CFIG4.eps

