basic_stats.f90

module basic_stats_mod
use iso_fortran_env, only: output_unit
use kind_mod, only: dp
use util_mod, only: default, print_table
use, intrinsic :: ieee_arithmetic, only: ieee_is_nan, ieee_value, ieee_quiet_nan
implicit none
private
public :: mean, variance, sd, mean_and_sd, kurtosis, basic_stats, &
   print_basic_stats, basic_stats_names, correl, acf, nbasic_stats, &
   stat, stats, corr_mat, rms, moving_sum, moving_average, moving_sd, moving_rms, &
   weighted_sd, &
   print_corr_mat, skew, cov, cov_mat, print_cov_mat, print_acf_mat, &
   print_acf, col_stats_ignore_nan, standardize_returns, biased_cov_sd
integer, parameter :: nbasic_stats = 6
character (len=*), parameter :: basic_stats_names(nbasic_stats) = &
   [character(len=4) :: "mean", "sd", "skew", "kurt", "min", "max"]
real(kind=dp), parameter :: bad_value = -huge(1.0d0)
interface stats
   module procedure stats_many_vec, stats_many_mat
end interface stats
interface print_acf
   module procedure print_acf_vec, print_acf_mat
end interface print_acf   
interface print_basic_stats
   module procedure print_basic_stats_vec, print_basic_stats_mat
end interface print_basic_stats
interface acf
   module procedure acf_vec, acf_mat
end interface acf
contains

pure function stats_many_vec(funcs, x) result(y)
! return statistics on x(:)
character (len=*), intent(in) :: funcs(:)
real(kind=dp), intent(in) :: x(:)
real(kind=dp) :: y(size(funcs))
integer :: i
do i=1,size(funcs)
   y(i) = stat(funcs(i), x)
end do
end function stats_many_vec

pure function stats_many_mat(funcs, x) result(y)
! return a matrix of statistics on each column of x(:,:)
character (len=*), intent(in) :: funcs(:)
real(kind=dp), intent(in) :: x(:,:)
real(kind=dp) :: y(size(funcs), size(x,2))
integer :: i
do i=1,size(x,2)
   y(:,i) = stats_many_vec(funcs, x(:,i))
end do
end function stats_many_mat

pure function stat(func, x) result(y)
! return a statistic on x(:)
character (len=*), intent(in) :: func
real(kind=dp), intent(in) :: x(:)
real(kind=dp)             :: y
select case(func)
   case ("mean")    ; y = mean(x)
   case ("sd")      ; y = sd(x)
   case ("variance"); y = variance(x)
   case ("skew")    ; y = skew(x)
   case ("kurt")    ; y = kurtosis(x)
   case ("min")     ; y = minval(x)
   case ("max")     ; y = maxval(x)
   case ("first")
      if (size(x) > 0) then
         y = x(1)
      else
         y = bad_value
      end if 
   case ("last")
      if (size(x) > 0) then
         y = x(size(x))
      else
         y = bad_value
      end if 
   case default ; y = -huge(x)
end select
end function stat

pure function mean(x) result(xmean)
! return the mean of x(:)
real(kind=dp), intent(in) :: x(:)
real(kind=dp) :: xmean
xmean = sum(x)/max(1,size(x))
end function mean

pure function sd(x) result(xsd)
! return the standard deviation of x(:)
real(kind=dp), intent(in) :: x(:)
real(kind=dp) :: xsd
real(kind=dp) :: m, var
integer :: n
n = size(x)
m = sum(x) / n
var = sum((x - m)**2) / (n-1)
xsd = sqrt(max(0.0_dp, var))
end function sd

pure function rms(x) result(xrms)
! return the root-mean-square of x(:)
real(kind=dp), intent(in) :: x(:)
real(kind=dp) :: xrms
xrms = sqrt(sum(x**2)/size(x))
end function rms

pure function mean_and_sd(x) result(res)
! return the mean and standard deviation of x(:)
real(kind=dp), intent(in) :: x(:)
real(kind=dp)             :: res(2)
real(kind=dp)             :: var
integer :: n
n = size(x)
res(1) = sum(x) / n
var = sum((x - res(1))**2) / (n-1)
res(2) = sqrt(max(0.0_dp, var))
end function mean_and_sd

pure function variance(x) result(var)
! return the variance of x(:)
real(kind=dp), intent(in) :: x(:)
real(kind=dp) :: var, m
integer :: n
n = size(x)
m = sum(x) / n
var = sum((x - m)**2) / (n-1)
end function variance

pure function skew(x) result(skew_val)
! return the skewness of x
real(kind=dp), intent(in) :: x(:)
real(kind=dp) :: skew_val
real(kind=dp) :: mean_x, sd_x
integer :: n
n = size(x)
mean_x = mean(x)
sd_x = sd(x)
skew_val = sum(((x - mean_x) / sd_x)**3) / n
end function skew

pure function kurtosis(x) result(kurtosis_val)
! return the kurtosis of x
real(kind=dp), intent(in) :: x(:)
real(kind=dp) :: kurtosis_val
real(kind=dp) :: mean_x, sd_x
integer :: n
n = size(x)
mean_x = mean(x)
sd_x = sd(x)
kurtosis_val = sum(((x - mean_x) / sd_x)**4) / n - 3.0_dp
end function kurtosis

pure function basic_stats(x) result(stats)
real(kind=dp), intent(in) :: x(:)
real(kind=dp)             :: stats(nbasic_stats)
stats = [mean(x), sd(x), skew(x), kurtosis(x), minval(x), maxval(x)]
end function basic_stats

subroutine print_basic_stats_vec(x, outu, fmt_header, fmt_trailer, &
   title, fmt_r, fmt_stats_names)
! print stats on a 1-D array
real(kind=dp), intent(in) :: x(:)
integer, intent(in), optional :: outu
character (len=*), intent(in), optional :: fmt_header, fmt_trailer, &
   title, fmt_r, fmt_stats_names
character (len=100) :: fmt_r_, fmt_stats_names_
integer :: i, outu_
if (present(fmt_stats_names)) then
   fmt_stats_names_ = fmt_stats_names
else
   fmt_stats_names_ = "(*(a10))"
end if
if (present(fmt_r)) then
   fmt_r_ = fmt_r
else
   fmt_r_ = "(*(f10.4))"
end if
outu_ = default(output_unit, outu)
if (present(fmt_header)) write (outu_, fmt_header)
if (present(title)) write (outu_, "(a)") title
write (outu_, fmt_stats_names_) (trim(basic_stats_names(i)), i=1,nbasic_stats)
write (outu_, fmt_r_) basic_stats(x)
if (present(fmt_trailer)) write (outu_, fmt_trailer)
end subroutine print_basic_stats_vec

subroutine print_basic_stats_mat(x, labels, outu, &
   fmt_header, fmt_trailer, title, fmt_cr, fmt_stats_names)
! print stats on a 2-D array
real(kind=dp), intent(in) :: x(:,:)
character (len=*), intent(in) :: labels(:)
integer, intent(in), optional :: outu
character (len=*), intent(in), optional :: fmt_header, fmt_trailer, &
   title, fmt_cr, fmt_stats_names
character (len=100) :: fmt_cr_, fmt_stats_names_
integer :: i, outu_
if (present(fmt_stats_names)) then
   fmt_stats_names_ = fmt_stats_names
else
   fmt_stats_names_ = "(*(a10))"
end if
if (present(fmt_cr)) then
   fmt_cr_ = fmt_cr
else
   fmt_cr_ = "(*(f10.4))"
end if
outu_ = default(output_unit, outu)
if (present(fmt_header)) write (outu_, fmt_header)
if (present(title)) write (outu_, "(a)") title
write (outu_, fmt_stats_names_) "", (trim(basic_stats_names(i)), i=1,nbasic_stats)
do i=1,size(x, 2)
   write (outu_, fmt_cr_) trim(labels(i)), basic_stats(x(:,i))
end do
if (present(fmt_trailer)) write (outu_, fmt_trailer)
end subroutine print_basic_stats_mat

pure function correl(x, y) result(corr_xy)
! Returns the linear Pearson correlation of x(:) and y(:)
! Returns a correlation < -1.0_dp to signal an error
real(kind=dp), intent(in) :: x(:), y(:)
real(kind=dp) :: corr_xy
real(kind=dp) :: x_mean, y_mean, cov_xy, var_x, var_y
integer :: n
n = size(x)
if (n /= size(y) .or. n == 0) then
   corr_xy = -2.0_dp
   return
end if
x_mean = sum(x) / n
y_mean = sum(y) / n
cov_xy = sum((x - x_mean) * (y - y_mean))
var_x  = sum((x - x_mean)**2)
var_y  = sum((y - y_mean)**2)
if (var_x <= 0.0_dp .or. var_y <= 0.0_dp) then
   corr_xy = -3.0_dp
else
   corr_xy = cov_xy / sqrt(var_x * var_y)
end if
end function correl

pure function cov(x, y) result(cov_xy)
! Returns the covariance of two 1D arrays
real(kind=dp), intent(in) :: x(:), y(:)
real(kind=dp) :: cov_xy
real(kind=dp) :: x_mean, y_mean
integer :: n
n = size(x)
if (n /= size(y) .or. n == 0) then
   error stop "x and y must have same size > 0 in cov"
end if
x_mean = sum(x) / n
y_mean = sum(y) / n
cov_xy = sum((x - x_mean) * (y - y_mean))
end function cov

pure function acf_vec(x, nacf) result(xacf)
! return the autocorrelations at lags 1 through nacf
real(kind=dp), intent(in) :: x(:)         ! Input array
integer, intent(in) :: nacf               ! Number of autocorrelations to compute
real(kind=dp) :: xacf(nacf)               ! Output array for autocorrelations
real(kind=dp) :: denom
real(kind=dp), allocatable :: xdm(:)      ! Demeaned version of x
integer :: n, lag
n = size(x)
xdm = x - mean(x)                          ! Compute demeaned x
denom = sum(xdm**2)
! Compute autocorrelation for each lag from 1 to nacf
do lag = 1, nacf
   xacf(lag) = sum(xdm(1:n-lag) * xdm(lag+1:n)) / denom
end do
end function acf_vec

pure function acf_mat(x, nacf) result(xacf)
! return the autocorrelations at lags 1 through nacf
real(kind=dp), intent(in) :: x(:,:)       ! Input array
integer, intent(in) :: nacf               ! Number of autocorrelations to compute
real(kind=dp) :: xacf(nacf,size(x,2))     ! Output array for autocorrelations
integer :: icol
do icol=1,size(x,2)
   xacf(:,icol) = acf(x(:,icol), nacf)
end do
end function acf_mat

subroutine print_acf_vec(x, nacf, label, outu, fmt_header, &
   fmt_trailer, title, fmt_acf, fmt_label)
! print the autocorrelations at lags 1 through nacf of x(:)
real(kind=dp), intent(in) :: x(:)       ! Input array
integer, intent(in) :: nacf             ! Number of autocorrelations to compute
character (len=*), intent(in), optional :: title, label, &
   fmt_header, fmt_trailer, fmt_acf, fmt_label
character (len=100) :: fmt_acf_, fmt_label_
integer, intent(in), optional :: outu
real(kind=dp) :: xacf(nacf)     
integer :: iacf, outu_
outu_ = default(output_unit, outu)
fmt_label_ = default("(6x,a8)", fmt_label)
if (present(fmt_header)) write (outu_, fmt_header)
if (present(title)) write (outu_, "(a)") title
if (present(label)) write (outu_,fmt_label_) label
fmt_acf_ = default("('ACF_', i2.2, f8.4)", fmt_acf)
xacf = acf_vec(x, nacf)
do iacf=1,nacf
   write (outu_, fmt_acf_) iacf, xacf(iacf)
end do
if (present(fmt_trailer)) write (outu_, fmt_trailer)
end subroutine print_acf_vec

subroutine print_acf_mat(x, nacf, labels, outu, fmt_header, &
   fmt_trailer, title, fmt_acf, fmt_labels)
! print the autocorrelations at lags 1 t hrough nacf of the columns of x(:,:)
real(kind=dp), intent(in) :: x(:,:)       ! Input array
integer, intent(in) :: nacf               ! Number of autocorrelations to compute
character (len=*), intent(in), optional :: title, labels(:), &
   fmt_header, fmt_trailer, fmt_acf, fmt_labels
integer, intent(in), optional :: outu
real(kind=dp) :: xacf(nacf,size(x,2))
integer :: iacf, icol, outu_
character (len=100) :: fmt_acf_, fmt_labels_
outu_ = default(output_unit, outu)
fmt_labels_ = default("(6x,*(a8))", fmt_labels)
if (present(fmt_header)) then
   write (outu_, fmt_header)
end if
if (present(title)) write (outu_, "(a)") title
if (present(labels)) write (outu_, fmt_labels_) &
   (trim(labels(icol)), icol=1,size(labels))
xacf = acf_mat(x, nacf)
fmt_acf_ = default("('ACF_', i2.2, *(f8.4))", fmt_acf)
do iacf=1,nacf
   write (outu_, fmt_acf_) iacf, xacf(iacf,:)
end do
if (present(fmt_trailer)) write (outu_, fmt_trailer)
end subroutine print_acf_mat

subroutine print_corr_mat(x, col_names, outu, fmt_col_names, fmt_row, &
   fmt_header, fmt_trailer)
! print the correlation matrix of the columns of x(:,:)
real(kind=dp), intent(in) :: x(:,:)
character (len=*), intent(in) :: col_names(:)
integer          , intent(in), optional :: outu ! output unit
character (len=*), intent(in), optional :: fmt_header, fmt_trailer, &
   fmt_col_names, fmt_row
character (len=100) :: fmt_col_names_, fmt_row_
fmt_col_names_ = default("(*(a8,:,1x))", fmt_col_names)
fmt_row_ = default("(a8, *(1x,f8.4))", fmt_row)
call print_table(corr_mat(x), row_names=col_names, col_names=col_names, &
   fmt_header=fmt_header, fmt_trailer=fmt_trailer, outu=outu, &
   fmt_col_names=fmt_col_names_, fmt_row=fmt_row_)
end subroutine print_corr_mat

subroutine print_cov_mat(x, col_names, outu, fmt_col_names, fmt_row, &
   fmt_header, fmt_trailer)
! print the covariance matrix of the columns of x(:,:)
real(kind=dp), intent(in) :: x(:,:)
character (len=*), intent(in) :: col_names(:)
integer          , intent(in), optional :: outu ! output unit
character (len=*), intent(in), optional :: fmt_header, fmt_trailer, &
   fmt_col_names, fmt_row
character (len=100) :: fmt_col_names_, fmt_row_
fmt_col_names_ = default("(*(a8,:,1x))", fmt_col_names)
fmt_row_ = default("(a8, *(1x,f8.4))", fmt_row)
call print_table(cov_mat(x), row_names=col_names, col_names=col_names, &
   fmt_header=fmt_header, fmt_trailer=fmt_trailer, outu=outu, &
   fmt_col_names=fmt_col_names_, fmt_row=fmt_row_)
end subroutine print_cov_mat

pure function corr_mat(x) result(cor)
    ! return the correlation matrix of the columns of x(:,:)
    real(kind=dp), intent(in) :: x(:,:)
    real(kind=dp)             :: cor(size(x,2), size(x,2))
    real(kind=dp)             :: mean_vec(size(x,2)), std_vec(size(x,2))
    real(kind=dp)             :: centered_x(size(x,1), size(x,2))
    integer                   :: n, p

    n = size(x, 1)  ! Number of rows
    p = size(x, 2)  ! Number of columns

    ! Compute the mean of each column
    mean_vec = sum(x, dim=1) / n

    ! Center the matrix by subtracting the mean of each column
    centered_x = x - spread(mean_vec, dim=1, ncopies=n)

    ! Compute the standard deviation of each column
    std_vec = sqrt(sum(centered_x**2, dim=1) / (n - 1))

    cor = matmul(transpose(centered_x), centered_x) / (n - 1)
    cor = cor / spread(std_vec, dim=1, ncopies=p)
    cor = cor / spread(std_vec, dim=2, ncopies=p)
end function corr_mat

pure function cov_mat(x) result(xcov)
    ! return the covariance matrix of the columns of x(:,:)
    real(kind=dp), intent(in) :: x(:,:)
    real(kind=dp)             :: xcov(size(x,2), size(x,2))
    real(kind=dp)             :: mean_vec(size(x,2)), std_vec(size(x,2))
    real(kind=dp)             :: centered_x(size(x,1), size(x,2))
    integer                   :: n, p

    n = size(x, 1)  ! Number of rows
    p = size(x, 2)  ! Number of columns

    ! Compute the mean of each column
    mean_vec = sum(x, dim=1) / n

    ! Center the matrix by subtracting the mean of each column
    centered_x = x - spread(mean_vec, dim=1, ncopies=n)

    ! Compute the standard deviation of each column
    std_vec = sqrt(sum(centered_x**2, dim=1) / (n - 1))
    xcov = matmul(transpose(centered_x), centered_x) / (n - 1)
end function cov_mat

pure function moving_sum(x, k) result(xsum)
! return a moving sum of x(:) with k terms, using fewer terms for i < k
real(kind=dp), intent(in) :: x(:)
integer      , intent(in) :: k
real(kind=dp)             :: xsum(size(x))
integer                   :: i, n
n = size(x)
if (n < 1) return
if (k < 1) then
   xsum = 0.0_dp
   return
end if
xsum(1) = x(1)
do i=2,min(k, n)
   xsum(i) = xsum(i-1) + x(i)
end do
do i=k+1, n
   xsum(i) = xsum(i-1) + x(i) - x(i-k)
end do
end function moving_sum

pure function moving_average(x, k) result(xma)
! return a moving average of x(:) with k terms, using fewer terms for i < k
real(kind=dp), intent(in) :: x(:)
integer      , intent(in) :: k
real(kind=dp)             :: xma(size(x))
integer                   :: i, n
real(kind=dp)             :: xsum(size(x))
n = size(x)
if (k < 1) then
   xma = 0.0_dp
   return
end if
xsum = moving_sum(x, k)
do i=1,min(k, n)
   xma(i) = xsum(i)/i
end do
do i=k+1,n
   xma(i) = xsum(i)/k
end do
end function moving_average

pure subroutine col_stats_ignore_nan(x, n, mean_x, sd_x, min_x, max_x) ! compute basic stats ignoring nan values
real(kind=dp), intent(in) :: x(:)
integer, intent(out) :: n
real(kind=dp), intent(out) :: mean_x, sd_x, min_x, max_x
real(kind=dp) :: s1, s2, xi, nan
integer :: i

nan = ieee_value(0.0_dp, ieee_quiet_nan)
n = 0
s1 = 0.0_dp
s2 = 0.0_dp
min_x = nan
max_x = nan

do i = 1, size(x)
   xi = x(i)
   if (ieee_is_nan(xi)) cycle
   n = n + 1
   s1 = s1 + xi
   s2 = s2 + xi*xi
   if (n == 1) then
      min_x = xi
      max_x = xi
   else
      if (xi < min_x) min_x = xi
      if (xi > max_x) max_x = xi
   end if
end do

if (n == 0) then
   mean_x = nan
   sd_x = nan
   min_x = nan
   max_x = nan
else if (n == 1) then
   mean_x = s1
   sd_x = nan
else
   mean_x = s1/n
   sd_x = sqrt(max(0.0_dp, (s2 - n*mean_x*mean_x)/(n - 1)))
end if
end subroutine col_stats_ignore_nan

pure function moving_sd(x, k) result(xsd)
! return a moving standard deviation of x(:) with k terms
real(kind=dp), intent(in) :: x(:)
integer      , intent(in) :: k
real(kind=dp)             :: xsd(size(x))
integer                   :: i, n
real(kind=dp)             :: xsum(size(x)), xsum2(size(x))
real(kind=dp)             :: m, v
n = size(x)
xsd = ieee_value(0.0_dp, ieee_quiet_nan)
if (k < 2 .or. n < k) return
xsum = moving_sum(x, k)
xsum2 = moving_sum(x**2, k)
do i=k, n
   m = xsum(i)/k
   v = (xsum2(i) - k*m**2)/(k-1)
   xsd(i) = sqrt(max(0.0_dp, v))
end do
end function moving_sd

pure function moving_rms(x, k) result(xrms)
! return a moving root-mean-square of x(:) with window k, assuming zero mean
real(kind=dp), intent(in) :: x(:)  ! input time series
integer      , intent(in) :: k     ! window length
real(kind=dp)             :: xrms(size(x))
integer                   :: i, n
real(kind=dp)             :: xsum2(size(x))
n = size(x)
xrms = ieee_value(0.0_dp, ieee_quiet_nan)
if (k < 1 .or. n < k) return
xsum2 = moving_sum(x**2, k)
do i = k, n
   xrms(i) = sqrt(xsum2(i)/k)
end do
end function moving_rms

pure function weighted_sd(x, w) result(xsd)
! return the weighted standard deviation of x(:) with weights w(:)
real(kind=dp), intent(in) :: x(:), w(:)
real(kind=dp)             :: xsd
real(kind=dp)             :: w_sum, xw_mean, var
integer                   :: n
n = size(x)
if (n < 2 .or. size(w) /= n) then
   xsd = ieee_value(0.0_dp, ieee_quiet_nan)
   return
end if
w_sum = sum(w)
if (w_sum <= 0.0_dp) then
   xsd = ieee_value(0.0_dp, ieee_quiet_nan)
   return
end if
xw_mean = sum(w * x) / w_sum
var = sum(w * (x - xw_mean)**2) / w_sum
xsd = sqrt(max(0.0_dp, var))
end function weighted_sd

function standardize_returns(R, use_ewma, ewma_lambda) result(R_std)
    !> Standardise a returns matrix R(n, p) column-by-column.
    !> use_ewma=.true.:  EWMA (RiskMetrics) normalisation with decay ewma_lambda.
    !> use_ewma=.false.: global normalisation (subtract mean, divide by full-sample sd).
    real(kind=dp), intent(in) :: R(:,:)
    logical, intent(in)       :: use_ewma
    real(kind=dp), intent(in) :: ewma_lambda
    real(kind=dp), allocatable :: R_std(:,:)
    real(kind=dp), allocatable :: glob_mean(:), glob_sd(:)
    real(kind=dp) :: var_t, sig_t
    integer :: n, p, a, i

    n = size(R, 1)
    p = size(R, 2)
    allocate(R_std(n, p), glob_mean(p), glob_sd(p))

    do a = 1, p
        glob_mean(a) = sum(R(:, a)) / n
        glob_sd(a)   = sqrt(sum((R(:, a) - glob_mean(a))**2) / n)
    end do

    if (use_ewma) then
        print "('normalisation: EWMA (lambda=',f0.2,')')", ewma_lambda
        do a = 1, p
            var_t = max(glob_sd(a)**2, tiny(1.0_dp))
            do i = 1, n
                sig_t      = sqrt(var_t)
                R_std(i,a) = (R(i,a) - glob_mean(a)) / sig_t
                var_t      = ewma_lambda * var_t + &
                             (1.0_dp - ewma_lambda) * (R(i,a) - glob_mean(a))**2
                var_t      = max(var_t, tiny(1.0_dp))
            end do
        end do
    else
        print "('normalisation: global (constant sigma)')"
        do a = 1, p
            if (glob_sd(a) > 0.0_dp) then
                R_std(:, a) = (R(:, a) - glob_mean(a)) / glob_sd(a)
            else
                R_std(:, a) = 0.0_dp
            end if
        end do
    end if
end function standardize_returns

pure subroutine biased_cov_sd(R, S, sd_vec)
    !> Compute the biased sample covariance matrix S and standard deviations sd_vec
    !! of the columns of R.  Divides by n (not n-1).
    real(kind=dp), intent(in)  :: R(:,:)
    real(kind=dp), intent(out) :: S(:,:), sd_vec(:)
    real(kind=dp) :: mu(size(R,2))
    integer :: a, b, m, p
    m = size(R, 1)
    p = size(R, 2)
    do a = 1, p
        mu(a) = sum(R(:, a)) / m
    end do
    do a = 1, p
        do b = a, p
            S(a,b) = sum((R(:,a) - mu(a)) * (R(:,b) - mu(b))) / m
            S(b,a) = S(a,b)
        end do
    end do
    do a = 1, p
        sd_vec(a) = sqrt(max(S(a,a), 0.0_dp))
    end do
end subroutine biased_cov_sd

end module basic_stats_mod