Option Base 1
Option Explicit

'matXは1元配列か2元配列かセル範囲のvariant型であって，その中にelementが含まれているなら，Trueを返す
Private Function isExistS(matX, element) As Boolean
With WorksheetFunction
    On Error GoTo Err:
        If Not .IsError(.Match(element, matX, 0)) Then
            isExistS = True
            Exit Function
        Else
            isExistS = False
            Exit Function
        End If
    On Error GoTo 0
    Exit Function
Err:
    isExistS = False
End With
End Function

'yを定数項とX0で重回帰．Excel Linest関数利用．Excelとは違い教科書通りに列ベクトルで返す
Function OLSE(y, X0)
With WorksheetFunction
    Dim bhatByExcel, k, i, bhat
    
    bhatByExcel = .LinEst(y, X0)
    k = UBound(bhatByExcel)
    '1元配列bhatByExcelを，列ベクトル化して逆順にする
    ReDim bhat(k, 1)
    For i = 0 To k - 1
        bhat(1 + i, 1) = bhatByExcel(k - i)
    Next
    OLSE = bhat
End With
End Function

'選択範囲で重複がなければTrue
Private Function isUniqueRng(rngX) As Boolean
    Dim element

    On Error GoTo Err
    For Each element In rngX
        If WorksheetFunction.CountIf(rngX, element) > 1 Then
            isUniqueRng = False
            Exit Function
        End If
    Next
    On Error GoTo 0
    
    isUniqueRng = True
    Exit Function
    
Err:
    isUniqueRng = False
End Function

'2019以前ではunique関数がないので作成．成分に重複がなければTrue
Private Function isUnique(ary) As Boolean
    Dim pivot, i, checkVal
    
    'Ubound(ary)にエラーが出るなら，その多くはaryがセル範囲のため．このときisUniqueRng関数を利用
    On Error GoTo Err0
        checkVal = UBound(ary)
    On Error Resume Next

    '以降はaryが1元配列
    On Error GoTo Err1
        For pivot = 1 To UBound(ary)
            For i = pivot + 1 To UBound(ary)
                If ary(i) = ary(pivot) Then
                    isUnique = False
                    Exit Function
                End If
            Next
        Next
    On Error GoTo 0
    
    isUnique = True
    Exit Function

Err0:
    isUnique = isUniqueRng(ary)
    Exit Function
Err1:
    isUnique = False
End Function

'aryの成分が全てminValue以上maxValue以下の整数なら，Trueを返す
Private Function isIntegerInRange(ary, minValue, maxValue) As Boolean
    Dim element
    
    On Error GoTo Err:
        For Each element In ary
            If element < minValue Or element > maxValue Or Int(element) <> element Then
                isIntegerInRange = False
                Exit Function
            End If
        Next
    On Error GoTo 0
    
    isIntegerInRange = True
    Exit Function
Err:
    isIntegerInRange = False
End Function

'F分布を使った変数を落とすF検定．第3引数以降に，第2引数の何列目を落とすかを書いていく．複数可
Function OLSdropVarsFtest(rngY, rngX0, ParamArray paramDropColumnNumberS())
With WorksheetFunction
    Dim y, X0, n, k0, k, J
    Dim olsResult0, olsResult1
    Dim aryDropColumnS, colX0, colX1
    Dim ssResidual0, ssResidual1, fNume, fDenom, fValue, pValue
    Dim strFvalue, strDf, strPvalue
    Dim i
    'y,n,k,JはOLSのものと同じ．X0は定数項のない説明変数行列．
    y = rngY.Value
    X0 = rngX0.Value
    n = UBound(X0, 1)
    k0 = UBound(X0, 2)
    k = k0 + 1
    J = UBound(paramDropColumnNumberS) + 1
    
    If J <= 0 Then
        OLSdropVarsFtest = "error:第3引数以降に，説明変数から落としたい変数の列番号の数字を入れてもう一度"
        Exit Function
    ElseIf J > k0 Then
        OLSdropVarsFtest = "error:第3引数以降に書いた変数の列番号の数字の個数は，全部でXの列数以下にしてもう一度"
        Exit Function
    End If
        
    '制約付きの回帰分析を行うために，落とす列の数字が入った配列を作る
    'dropColumnNumberSの配列は常に0スタートなのでこの表現となる．
    ReDim aryDropColumnS(J)
    For i = 0 To J - 1
        aryDropColumnS(1 + i) = paramDropColumnNumberS(i)
    Next
    'エラーチェック
    If Not isIntegerInRange(aryDropColumnS, 1, k0) Then
        OLSdropVarsFtest = "error:第3引数以降には，全て説明変数から落としたい変数の列番号の正整数を入れてもう一度"
        Exit Function
    End If
    If Not isUnique(aryDropColumnS) Then
        OLSdropVarsFtest = "error:第3引数以降に，重複があります．変更してもう一度"
        Exit Function
    End If
        
    '制約なしの回帰分析の残差
    olsResult0 = OLSEst3aryS(y, X0)
    ssResidual0 = .SumSq(olsResult0(3))
        
    '制約ありの回帰分析の残差
    'Jがk0の場合は，残差は総変動で，J<k0のときは制約付きの回帰から残差を作成
    If J = k0 Then
        ssResidual1 = .DevSq(y)
    Else
        'X0から，指定された列を落としたX1を作る
        ReDim X1(n, k0 - J)
        colX1 = 1
        For colX0 = 1 To k0
            If Not isExistS(aryDropColumnS, colX0) Then
                For i = 1 To n
                    X1(i, colX1) = X0(i, colX0)
                Next
                colX1 = colX1 + 1
            End If
        Next
        olsResult1 = OLSEst3aryS(y, X1)
        ssResidual1 = .SumSq(olsResult1(3))
    End If
        
    '2つの残差をもとに，F値とP値を作成
    fNume = (ssResidual1 - ssResidual0) / J
    fDenom = ssResidual0 / (n - k)
    fValue = fNume / fDenom
    pValue = .F_Dist_RT(fValue, J, n - k)
    
    '文字列を作成して，返す
    strFvalue = "F値=" & .Round(fValue, 4) & "，"
    strDf = "自由度(" & J & "," & n - k & ")のF分布より"
    strPvalue = "P値=" & .Round(pValue, 4)
    
    OLSdropVarsFtest = strFvalue & strDf & strPvalue
End With
End Function

'OLS推定．bhat,yhat,ehatの3つを返す．
Private Function OLSEst3aryS(y, X0)
With WorksheetFunction
    Dim n, k, bhat, X, yhat, ehat
    Dim i, jCol
    'n,k,bhatは，OLSと同じ
    bhat = OLSE(y, X0)
    n = UBound(X0, 1)
    k = UBound(X0, 2) + 1
    '定数項付きのXを作成
    ReDim X(n, k)
    For i = 1 To n
        X(i, 1) = 1
        For jCol = 2 To k
            X(i, jCol) = X0(i, jCol - 1)
        Next
    Next
    'OLSのyhatとehatを作成
    yhat = .MMult(X, bhat)
    ReDim ehat(n, 1)
    For i = 1 To n
        ehat(i, 1) = y(i, 1) - yhat(i, 1)
    Next
    
    OLSEst3aryS = Array(bhat, yhat, ehat)
End With
End Function
