Merge pull request #23 from atecon/22-wrong-p-value-for-cg-sign-test

correct pvalue for the center of CGSIGN test dist

docs/FEP.lyx

@@ -131,7 +131,7 @@ based on FEP
 
 \end_inset
 
- version 2.10, April 2023
+ version 2.91, May 2023
 \end_layout
 
 \begin_layout Abstract
@@ -14150,7 +14150,7 @@ Changelog
 \end_layout
 
 \begin_layout Itemize
-version 2.10 (May 2023)
+version 2 .91 (May 2023)
 \end_layout
 
 \begin_deeper
@@ -14167,6 +14167,14 @@ extra
  version is 0.41)
 \end_layout
 
+\begin_layout Itemize
+Bugfix: Computing two-sided Sign test for 
+\family typewriter
+doCGtest()
+\family default
+ was wrong in some case
+\end_layout
+
 \begin_layout Itemize
 Minor internal refactorings
 \end_layout

run_tests.sh

@@ -4,16 +4,18 @@ set -e
 DIR=$(dirname $(realpath "$0")) 	# locate folder where this sh-script is located in
 
 PROJECT="FEP"
-SCRIPT_1="./tests/test_forecast_metrics.inp"
+
+FILELIST="./tests/test_forecast_metrics.inp \
+          ./tests/test_private.inp"
 
 
 cd $DIR
 echo "Switched to ${DIR}"
 
 
-# Execute scripts
-gretlcli -b -e -q ${SCRIPT_1}
-
+for f in ${FILELIST};
+  do gretlcli -b -e -q ${f}
+done;
 
 if [ $? -eq 0 ]
 then

src/FEP.spec

@@ -1,7 +1,7 @@
 author = Artur Tarassow and Sven Schreiber
 email = atecon@posteo.de
 version = 2.91
-date = 2023-05-03
+date = 2023-05-04
 description = Forecast Evaluation Package
 tags = C12 C52 C53
 min-version = 2021a

src/FEP_private.inp

@@ -479,6 +479,55 @@ function matrix rattest (matrix Y,
 end function
 
 
+function scalar pvalue_two_sided_signtest (const int T[1::],
+                                           const scalar teststat)
+    /* Computes the p-value for a sign test using a two-sided alternative.
+
+    Args:
+        T (int]): the number of trials.
+        teststat (scalar): the sign test statistic.
+
+    The forecast error is E = Y - F.
+    Under the null H0: p=0.5 (equal probability which implies that the median of E is zero.
+
+    W is the number of pairs that E > 0. If H0 is true, then W follows a
+    binomial distribution W ~ b(T, 0.5).
+
+    The left-tail value is computed by Pr(W ≤ w), which is the p-value for the
+    alternative H1: p < 0.50. This alternative means that the F measurements
+    tend to be higher.
+
+    The right-tail value is computed by Pr(W ≥ w), which is the p-value for
+    the alternative H1: p > 0.50. This alternative means that the Y
+    measurements tend to be higher.
+
+    For a two-sided alternative H1 the p-value is twice the smaller tail-value.
+
+    The C&G test is two-sided: H0: Median of E = 0, H1: Median of E != 0
+
+    Source: https://en.wikipedia.org/wiki/Sign_test
+    */
+
+    if teststat > T
+        printf "Warning: stat %d out of range for %d trials!\n", teststat, T
+        return NA
+    elif teststat < T/2
+        # left tail of the distribution (for odd and even T)
+        # (gretl does X <= x with cdf())
+        # because of symmetry under p=0.5
+        scalar pval = 2 * cdf(B, 0.5, T, teststat)
+    elif teststat == T/2    # center
+        scalar pval = 1
+    else                    # right tail
+        # (force X >= x by passing x' = x-1, since gretl does
+        # X > x with pvalue) because of symmetry under p=0.5
+        scalar pval = 2 * pvalue(B, 0.5, T, teststat - 1)
+    endif
+
+    return pval
+end function
+
+
 /*
 # not needed anymore, replaced by native lrcovar
 

src/FEP_public.inp

@@ -394,9 +394,9 @@ end function
 
 # ---------------------------------------------
 
-function matrices CamDufStats( matrix mdat,
-                              matrix whichlags[null], int recenter[0:2:0])
-    # (switched the arg ordering ...)
+function matrices CamDufStats (matrix mdat,
+                               matrix whichlags[null],
+                               int recenter[0:2:0])
     /*
        Function to calculate the test statistics from
        Campbell & Dufour 1995 (C&D), which are sign statistics

src/FEP_tests_individual.inp

@@ -499,17 +499,7 @@ function void doCGtest (bundle *b)
     ## sign test
     scalar b.CGSIGNstat = testout[1][1,1]
 
-    # Binomial dist
-    if b.CGSIGNstat < (T + 1) / 2 # left tail of the distribution
-        temp = pvalue(B, 0.5, T, b.CGSIGNstat)
-        # (gretl does > x, not >= x here!)
-        b.CGSIGNpval = 2 * (1 - temp) # because of symmetry under p=0.5
-
-    else	# right tail
-        # (force X >= x by passing x' = x-1, since gretl does > x)
-        b.CGSIGNpval = 2 * pvalue(B, 0.5, T, b.CGSIGNstat - 1)
-        # (again, because of symmetry)
-    endif
+    b.CGSIGNpval = pvalue_two_sided_signtest(T, b.CGSIGNstat)
 
     # signed rank test
     if !inbundle(b, "CGX") # serial correlation, or unbiasedness

tests/test_private.inp

@@ -0,0 +1,40 @@
+set verbose off
+clear
+
+set assert stop
+include assertion.gfn
+
+include "./src/FEP.inp" --force
+
+
+bundles Params = null
+Params = Params + _(T=100, teststat = 50, expected = 1)
+Params = Params + _(T=100, teststat = 101, expected = NA)
+Params = Params + _(T=100, teststat = 25,
+                    expected = 2 * cdf(B, 0.5, 100, 25))
+Params = Params + _(T=100, teststat = 60,
+                    expected = 2 * pvalue(B, 0.5, 100, 60 - 1))
+
+Params = Params + _(T=10, teststat = 8, expected = 0.109375)
+
+function void test_pvalue_two_sided_signtest (const bundles P)
+    print "Start testing function pvalue_two_sided_signtest()."
+
+    loop foreach i P
+        # Given + When
+        scalar actual = pvalue_two_sided_signtest(P[i].T, P[i].teststat)
+
+        # Then
+        if ok(actual) && ok(P[i].expected)
+            assert_almost_equal_num(actual, P[i].expected, 10e-7)
+        elif !ok(actual) && !ok(P[i].expected)
+            print "assertion is OK"
+        else
+            assert(TRUE == FALSE)
+        endif
+    endloop
+end function
+test_pvalue_two_sided_signtest(Params)
+
+
+printf "\nInfo: All tests passed.\n"