% [moments,reg,phiv]=NIGLPmoments(alpha,beta,delta,r,t,
%                   derivwanted, apscalwanted, plotwanted)
%
% This function calculates the µ-centered r-th moments
% of an NIG(alpha,beta,mu,delta) Levy-process at the times t 
% and, if wanted, the derivatives of the log µ-centered r-th 
% moment with respect to the log of time and the 
% apparent scaling coefficient. 
%
% Input:
%
% alpha: 
% parameter alpha of the NIG Levy process, positive 
% beta: 
% parameter beta of the NIG Levy process, less than alpha in absolute value
% delta: 
% parameter delta of the NIG Levy process, positive
% r: 
% order of the µ-centered  moment to be calculated, positive integer
% t: 
% time(s) at which the moments (and derivatives) are to be calculated, row vector
% derivwanted: 
% Boolean variable indicating whether or not derivatives are to be calculated (0= no derivatives)
% apscalwanted: 
% Boolean variable indicating whether or not the apparent scaling coefficient is to be calculated via regression (0= no)
% plotwanted: 
% Boolean variable indicating whether or not moments and regression line are to be plotted (0= no plots)
%
% Output:
%
% moments: 
% row vector of the moments (at times t)
% reg: 
% column vector of regression results (second entry: apparent scaling slope, first: intercept)
% phiv: 
% vector of values of the derivative of the log  µ-centered moment with respect to log time (at times t)
%
% For details on the formulae implemented see the paper 
% "Absolute moments of Generalized Hyperbolic Distributions 
% and Approximate Scaling of NIG Levy processes" 
% by Ole Eiler Barndorff-Nielsen (University of Aarhus) 
% and Robert Stelzer (Munich University of Technology) 
% 
% This code has been implemented in Matlab by Robert Stelzer. 
% It does not require any special toolboxes. It may be freely 
% used and distributed free of charge.
% It may be altered to fit different 
% needs upon the condition that there remains a reference 
% to the original author and the revised code is made available 
% free of charge.
%
% The author would be pleased to be informed about any research 
% results obtained using this code. (email: rstelzer (at) ma (dot) tum (dot) de)
%
% Robert Stelzer, 18.4.2004

function [moments,regcoef,phiv]=NIGLPmoments(alpha,beta,delta,r,t,derivwanted,apscalwanted, plotwanted)
%initializing
regcoef=0;
phiv=0;
n=length(t);  
rbar=ceil(r/2)-1/2;
remain=mod(r,2);
moments=zeros(1,n);
if (derivwanted)
phiv=zeros(1,n);
end;
% calculate moments 
for s=1:n
    %setup parameter gamma
    gam=sqrt(alpha*alpha-beta*beta);
    %setup bar parameters needed in the calculations
    alphabar=alpha*delta;
    gammabar=gam*delta;
    betabar=beta*delta;  
    talphabar=t(s)*alphabar;
    talphabarinv=1./talphabar;
    %initialize recursions
    m=2*betabar*betabar*t(s)/alphabar;
    a=(2*delta*delta*t(s)/alphabar).^(ceil(0.5*r)).*sqrt(2*t(s)*alphabar).*exp(t(s)*gammabar).*gamma(rbar+1)*beta^remain/pi;
    K0=besselk(rbar-1,talphabar);
    K1=besselk(rbar,talphabar);
    moments(s)=a.*K1;
    %initialize variables controlling recursions
    cont=1;
    k=1;
%compute coefficients and add them up, until neglegible (indicated by cont)
while cont
    a=(a.*m)*(k+rbar)/((2*k-1+remain)*(2*k+remain));
    Kt=K1;
    if K1==0 
        K1=besselk(k+rbar,talphabar);
    else
        K1=2*(k-1+rbar)*(talphabarinv.*K1)+K0;
    end;
    K0=Kt;
    add=a.*K1;
    if moments(s)~=0 
        cont=(add./moments(s)>1e-12);
    end;
    moments(s)=moments(s)+add;
    k=k+1;
end;
%Calculate dE|Z(t)-µt|/dln t if wanted
if (derivwanted)
  a=(2*delta*delta*t(s)/alphabar).^(ceil(0.5*r)).*sqrt(2*t(s)*alphabar).*exp(t(s)*gammabar).*gamma(rbar+1)*beta^remain/pi;
  K0=besselk(rbar-1,talphabar);
  K1=besselk(rbar,talphabar);
  numerator=a.*K0;
  cont=1;
  k=1;
    %compute coefficients and add them up, until neglegible (indicated by cont)
    while cont
        a=(a.*m)*(k+rbar)/((2*k-1+remain)*(2*k+remain));
        Kt=K1;
        if K1==0
            K1=besselk(k+rbar,talphabar)
        else
            K1=2*(k-1+rbar)*(talphabarinv.*K1)+K0;
        end;
        add=a.*K0;
        K0=Kt;
        if numerator~=0
            cont=(add/numerator>1e-12);
        end;
        numerator=numerator+add ;
        k=k+1;
    end;
    numerator
    moments(s)
    phiv(s)=1+gammabar*t(s)-alphabar*t(s).*numerator./moments(s);
end;
end;
if (apscalwanted)
    %regress logmoments vs logtime and constant
    logtime=zeros(length(t),1);
    logtime(:,1)=log(t)';
    logmoments=log(abs(moments));
    regmat=ones(length(t),2);
    regmat(:,2)=logtime;
    regcoef=regmat\logmoments';
    if (plotwanted)
        %plot log moments against log time together with regression matrix
         plot(logtime,logmoments,'.',logtime,regcoef(1)+regcoef(2)*logtime);
         xlabel('ln(t) in seconds');
         ylabel(strcat('ln ( |E(Z(t)-t\mu)^',sprintf('{%1.1f}',r),'|)'));
     end;
 end;