' ******************************************************************
' Test of weighted regression
' ******************************************************************

' ------------------------------------------------------------------
' Machine-dependent constant
' ------------------------------------------------------------------

#DEFINE MaxNum 1.797693134862315D+308
' Max. floating point number: 2^1024

' ------------------------------------------------------------------
' Type definition
' ------------------------------------------------------------------

TYPE TRegTest
  Vr   AS DOUBLE   ' Residual variance
  R2   AS DOUBLE   ' Coefficient of determination
  R2a  AS DOUBLE   ' Adjusted coeff. of determination
  F    AS DOUBLE   ' Variance ratio (explained/residual)
  Nu1  AS INTEGER  ' Number of d.o.f. (numerator)
  Nu2  AS INTEGER  ' Number of d.o.f. (denominator)
END TYPE

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from last function call

' ******************************************************************

SUB WRegTest(      Y()     AS DOUBLE, _
                   Ycalc() AS DOUBLE, _
                   S()     AS DOUBLE, _
                   V()     AS DOUBLE, _
             BYREF Test    AS TRegTest)
' ------------------------------------------------------------------
' Test of weighted regression
' ------------------------------------------------------------------
' Input parameters : Y()     = observed values of dependent variable
'                    Ycalc() = predicted values
'                    S()     = standard deviations
'                    V()     = inverse matrix
' Output parameters: V()     = variance-covariance matrix
'                    Test    = test results
' ------------------------------------------------------------------

  DIM AS INTEGER LbY, UbY, LbV, UbV, I, J, K, Nobs, Npar
  DIM AS DOUBLE  D, SSt, SSe, SSr, SW, SWY, Ybar

  LbY = LBOUND(Y)
  UbY = UBOUND(Y)

  LbV = LBOUND(V)
  UbV = UBOUND(V)

  Nobs = UbY - LbY + 1
  Npar = UbV - LbV + 1

  IF Nobs <= Npar THEN ErrCode = -1 : EXIT SUB

  DIM AS DOUBLE W(LbY TO UbY)  ' Weights
  
  SW = 0
  SWY = 0

  FOR K = LbY TO UbY
    IF S(K) <= 0 THEN ErrCode = -2 : EXIT SUB
    W(K) = 1 / (S(K) * S(K))
    SW = SW + W(K)
    SWY = SWY + W(K) * Y(K)
  NEXT K

  ErrCode = 0

  Ybar = SWY / SW

  SSt = 0
  SSe = 0
  SSr = 0

  FOR K = LbY TO UbY
    D = Y(K) - Ybar
    SSt = SSt + W(K) * D * D
    D = Ycalc(K) - Ybar
    SSe = SSe + W(K) * D * D
    D = Y(K) - Ycalc(K)
    SSr = SSr + W(K) * D * D
  NEXT K

  Test.Nu1 = Npar - 1
  Test.Nu2 = Nobs - Npar
  Test.R2 = SSe / SSt
  Test.R2a = 1 - (1 - Test.R2) * (Nobs - 1) / Test.Nu2
  Test.Vr = SSr / Test.Nu2

  IF SSr = 0 THEN
    Test.F = MaxNum
    EXIT SUB
  END IF

  Test.F = SSe / Test.Nu1 / Test.Vr

  ' Update var-cov matrix
  FOR I = LbV TO UbV
    FOR J = I TO UbV
      V(I, J) = V(I, J) * Test.Vr
      IF I <> J THEN V(J, I) = V(I, J)
    NEXT J
  NEXT I
END SUB

