Kod:
' Copyright 2002 Finnish Meteorological Institute
' Timo Salmi & Anu Määttä & Toni Amnell
' MAKESENS Version 1.0
'COMMON SETTINGS
Option Base 1 'Default indexing of arrays starts from 1
'Code for missing data in arrays:
Const MissingValue As Double = -999999#
'Maximum number of data in one time series:
Const MaxData As Integer = 100
'Codes for different significance levels:
'Minimum count of data to use normal approximation
' in Mann-Kendall test. Below this value the S statistics
' is used:
Const MinMannKendNorm As Integer = 10
'Minimum count of data to calculate confidence interval
'in Sen's method
Const MinSenConf As Integer = 10
Const S001 As String = "***" 'alpha = 0.001
Const S01 As String = "**" 'alpha = 0.01
Const S05 As String = "*" 'alpha = 0.05
Const S1 As String = "+" 'alpha = 0.1
'Arrays of critical values of Mann-Kendall statistic S
' for significance levels 0.001, 0.01, 0.05 and 0.1
' of two-sided test when n is between 4 and 10.
' The arrays will be filled by the subroutine fillS
Dim S_001(4 To 10) As Integer
Dim S_01(4 To 10) As Integer
Dim S_05(4 To 10) As Integer
Dim S_1(4 To 10) As Integer
Private Sub CB_CalculateStatistics_Click()
' The main program of calculation
' - Retrieves the data values from the sheet "Annual data" with
' the subroutine GetData
' - Calculates the statistics with the subroutines MannKendall
' and Sen and with the function calcB and saves the results into
' the sheet "Trend statistics"
' - Finally calls the workbook level routines makeCollection and
' drawFigure preparing the sheet Figure
Dim nofCol As Integer 'Number of columns i.e. time series
Dim colno As Integer 'Column number of a time series
Dim firstYear As Integer 'first year of a time series
Dim baseYear As Integer 'first year of all time series
Dim nYears As Integer 'number of years in a time series
Dim n As Integer 'true data values in a time series
'i.e. missing values are not considered
Dim x(MaxData) As Double 'Array for data values of a time series
Dim s As Integer 'Mann-Kendall test statistic for n=4..10
Dim Z As Double 'Mann-Kendall test statistic for n>10
Dim signif As String 'significance of trend
'Sen's slope estimator Q and its 99% and 95 % confidence levels:
Dim Q As Double, Qmin99 As Double, Qmax99 As Double
Dim Qmin95 As Double, Qmax95 As Double
'Constants B for equation of lines of Sen's slope and conf. intervals:
Dim B As Double, Bmin99 As Double, Bmax99 As Double
Dim Bmin95 As Double, Bmax95 As Double
' The result cells are emptied before the calculation starts
Worksheets("Trend Statistics").Range("E6:Q30") = ""
nofCol = Worksheets("Annual data").Cells(8, 2).value
baseYear = Worksheets("Annual data").Cells(14, 1).value
Call fillS 'initializes arrays of Mann-Kendall propabilities
'Calculation of trend statistics for each time series at a time
For colno = 2 To nofCol + 1
If Not GetData(colno, baseYear, firstYear, nYears, n, x) Then
Exit Sub
End If
If n >= 2 Then 'nothing can be computed, if n<2
'First the existence of trend is tested using Mann-Kendall method.
Call MannKendall(nYears, x, s, Z, signif)
If n < MinMannKendNorm Then
Worksheets("Trend Statistics").Cells(4 + colno, 5) = s
Else
Worksheets("Trend Statistics").Cells(4 + colno, 6) = Z
End If
Worksheets("Trend Statistics").Cells(4 + colno, 7) = signif
'Evaluation of Sen's slope estimator and confidence intervals
Call Sen(nYears, x, Q, Qmin99, Qmax99, Qmin95, Qmax95)
Worksheets("Trend Statistics").Cells(4 + colno, 8) = Q
B = calcB(nYears, x, firstYear, baseYear, Q)
Worksheets("Trend Statistics").Cells(4 + colno, 13) = B
If n >= MinSenConf Then
Worksheets("Trend Statistics").Cells(4 + colno, 9) = Qmin99
Worksheets("Trend Statistics").Cells(4 + colno, 10) = Qmax99
Worksheets("Trend Statistics").Cells(4 + colno, 11) = Qmin95
Worksheets("Trend Statistics").Cells(4 + colno, 12) = Qmax95
'Coefficients B for equation of linear trend f(t)=Qt+B
Bmin99 = calcB(nYears, x, firstYear, baseYear, Qmin99)
Bmax99 = calcB(nYears, x, firstYear, baseYear, Qmax99)
Bmin95 = calcB(nYears, x, firstYear, baseYear, Qmin95)
Bmax95 = calcB(nYears, x, firstYear, baseYear, Qmax95)
Worksheets("Trend Statistics").Cells(4 + colno, 14) = Bmin99
Worksheets("Trend Statistics").Cells(4 + colno, 15) = Bmax99
Worksheets("Trend Statistics").Cells(4 + colno, 16) = Bmin95
Worksheets("Trend Statistics").Cells(4 + colno, 17) = Bmax95
End If
End If
Next colno
' Draw the figure of the first component
Sheets("Figure").Cells(9, 3).value = 1
Sheets("Figure").Cells(10, 3).value = Sheets("Annual data").Cells(13, 2).value
Application.Run "makeCollection"
Application.Run "DrawFigure"
End Sub 'CB_CalculateStatistics_Click
Private Function GetData(ByVal colno As Integer, ByVal baseYear As Integer, firstYear As Integer, nYears As Integer, n As Integer, x() As Double) As Boolean
' Retrieving of data of one time series into the array x()
' colno is the column of the worksheet "Annual data" where the
' values of the time series exist.
' The real number of annual values n in time series is calculated
' If the cell is empty it is understood as a missing value.
Dim rowno As Integer 'row of the data cell
Dim lastYear As Integer 'last year of the time series
Dim nVal As Integer 'counter for number of true data
Dim i As Integer 'counter for data loop
Dim Error As Integer
firstYear = Worksheets("Annual data").Cells(10, colno).value
lastYear = Worksheets("Annual data").Cells(11, colno).value
nYears = lastYear - firstYear + 1
nVal = 0
For i = 1 To nYears
If firstYear < baseYear Then
Error = MsgBox("For the time series """ + _
Worksheets("Annual data").Cells(13, colno).value + _
""" first year is too early!")
GetData = False
Exit Function
End If
rowno = 13 + i + firstYear - baseYear
If IsEmpty(Worksheets("Annual data").Cells(rowno, colno)) Then
x(i) = MissingValue
Else
nVal = nVal + 1
x(i) = Worksheets("Annual data").Cells(rowno, colno)
End If
Next i
n = nVal
GetData = True
End Function ' GetData
Private Sub MannKendall(ByVal nYears As Integer, x() As Double, s As Integer, Z As Double, signif As String)
'Calculates the MannKendall test
'Calls the function tiedSum
'Uses the string constants S001, S01, S05 and S1
Dim absS As Integer 'value of absS
Dim varS As Double 'the variance of S
Dim absZ As Double 'value of abs(Z)
Dim k As Integer, j As Integer 'counters for slopes
Dim n As Long 'number of true values in x()
Z = MissingValue ' returns MissingValue for Z
' if they are not calculated
'Computing of the Mann-Kendall statistic S.
signif = ""
n = IIf(x(nYears) <> MissingValue, 1, 0)
s = 0
For k = 1 To nYears - 1
If x(k) <> MissingValue Then
n = n + 1
For j = k + 1 To nYears
If x(j) <> MissingValue Then
s = s + Sgn(x(j) - x(k))
End If
Next j
End If
Next k
If n < 4 Then
'If n is less than 4, the method can not be used at all
Exit Sub
ElseIf n < MinMannKendNorm Then
'If n is between 4 and 10, S is compared directly to Mann-Kendall statistics for S
absS = Abs(s)
signif = Switch(absS >= S_001(n), S001, absS >= S_01(n), S01, absS >= S_05(n), S05, absS >= S_1(n), S1, True, "")
Else 'n>=MinMannKendNorm
'If n is at least 10, the normal distribution is used
'Firstly the variance VAR(S) is calculated
'The correction term for ties is calculated by the function tiedSum
varS = (n * (n - 1) * (2 * n + 5) - tiedSum(nYears, x)) / 18#
'Calculation of test statistic Z using S and its variance VAR(S)
Z = Switch(s > 0, (s - 1) / Sqr(varS), s < 0, (s + 1) / Sqr(varS), s = 0, 0#)
'The absolute value of Z is compared to critical value Z[1-alpha/2]
'which is obtained from the standard normal table. The presence and
'significance of the trend is evaluated by testing four different
'levels of significance: '0.001, 0.01, 0.05 and 0.1
absZ = Abs(Z)
signif = Switch(absZ > 3.292, S001, absZ > 2.576, S01, absZ > 1.96, S05, absZ > 1.645, S1, True, "")
End If
End Sub 'MannKendall
Private Sub Sen(ByVal nYears As Integer, x() As Double, Q As Double, Qmin99 As Double, Qmax99 As Double, Qmin95 As Double, Qmax95 As Double)
'Calculates Sen's slope estimator Q and its 99% (Qmax99,Qmin99)
' and 95 % (Qmax95, Qmin95)confidence levels
' Calls the function tiedSum and the subroutine CalulateConfidenceInterval
Dim nofQ As Integer 'number of value pairs
Dim Qarray(MaxData * (MaxData - 1) / 2) As Double 'Array for the slopes of value pairs
Dim k As Integer, j As Integer 'counters for loops
Dim n As Long 'number of true values in x()
Dim Calpha As Double 'C-alpha for calculation of conf.intervals of Q
'Computing of slopes of individual value pairs into Qarray
nofQ = 0 'used as counter for Qarray
n = IIf(x(nYears) = MissingValue, 0, 1)
For k = 1 To nYears - 1
If x(k) <> MissingValue Then
n = n + 1
For j = k + 1 To nYears
If x(j) <> MissingValue Then
nofQ = nofQ + 1
Qarray(nofQ) = (x(j) - x(k)) / (j - k)
End If
Next j
End If
Next k
'The median of individual slopes in Qarray is the Sen's
'slope estimator. The median is calculated by the function "median".
Q = median(nofQ, Qarray)
If n >= MinSenConf Then
'The confidence intervals are calculated only if n is at least 10.
'Computing of variance VAR(S) of Mann-Kendall statistics S.
'The correction term for ties is calculated by the function tiedSum
varS = (n * (n - 1) * (2 * n + 5) - tiedSum(nYears, x)) / 18#
'The 100(1-alpha)% two-sided confidence intervals for the
'Sen's slope are computed with two values of alpha: 0.01 and 0.05
'which means 99 % and 95 % confidence intervals. The values of
'Z[1-alpha/2] are obtained from the standard normal table.
'Case alpha=0.01: Z[1-alpha/2]=Z[0.995]=2.576
Calpha = 2.576 * Sqr(varS)
Call CalculateConfidenceInterval(Calpha, nofQ, Qarray, Qmin99, Qmax99)
'Case alpha=0.05: Z[1-alpha/2]=1.96
Calpha = 1.96 * Sqr(varS)
Call CalculateConfidenceInterval(Calpha, nofQ, Qarray, Qmin95, Qmax95)
Else
Qmin99 = MissingValue
Qmax99 = MissingValue
Qmin95 = MissingValue
Qmax95 = MissingValue
End If
End Sub 'Sen
Private Function tiedSum(n As Integer, x() As Double) As Integer
'Calculates sum related to tied groups(= two or more equal values)
' for the variance of Mann-Kendall statistics S
'n = number of values in the array x including missing values
'Function tiedSum is called by subroutines Sen and MannKendallNorm
Dim m As Integer ' number of tied groups
Dim tval() As Double ' data values of tied groups
ReDim tval(n)
Dim t() As Integer, nt As Integer ' number of data in tied groups
ReDim t(n)
Dim p, i As Integer 'indexes for the loops
Dim newValue As Boolean
Dim tSum As Integer
'Calculation of the number of tied groups m and the number of data
' in tied groups t()
m = 0
For i = 1 To n - 1
If x(i) <> MissingValue Then
newValue = True
If m > 0 Then
For p = 1 To m
If x(i) = tval(p) Then
newValue = False 'this value is alredy managed
Exit For
End If
Next p
End If
If newValue Then
nt = 1 'number of equal values x(i)
For p = i + 1 To n
If x(p) = x(i) Then
nt = nt + 1
End If
Next p
If nt > 1 Then ' new group only if nt>1
m = m + 1
t(m) = nt
tval(m) = x(i)
End If
End If
End If
Next i
'Calculating the sum related to tied groups for variance
tSum = 0
If m > 0 Then
For p = 1 To m
tSum = tSum + t(p) * (t(p) - 1) * (2 * t(p) + 5)
Next p
End If
tiedSum = tSum
End Function 'tiedSum
Sub CalculateConfidenceInterval(ByVal Calpha As Double, ByVal nofQ As Integer, Qarray() As Double, lowerLimit As Double, upperLimit As Double)
'Computes confidence interval for Sen's slope estimate.
'Input parameters: Calpha = Z[1-alpha/2],
' nofQ - number of slopes of all data pairs
' Qarray - array of slopes of all data pairs
'Subroutine returns the lowerLimit and upperLimit.
'Calls the subroutine SortArray
'Is called by the suroutine Sen
Dim M1 As Double 'M1:th largest ordered slope
Dim M2 As Double 'M2:th largest ordered slope
Dim M1int As Integer 'integer part of M1 (>0)
Dim M2int As Integer 'integer part of M2+1 (>0)
Dim QarraySort() As Double
ReDim QarraySort(nofQ)
'The array Qarray is sorted to the array QarraySort
Call SortArray(nofQ, Qarray, QarraySort)
M1 = (nofQ - Calpha) / 2
M2 = (nofQ + Calpha) / 2
If M1 > 1 Then
'to be sure that index does not point outside QarraySort
M1int = Int(M1) 'find the integer part of M1
'Interpolation of the lower limit
lowerLimit = QarraySort(M1int) + (M1 - M1int) * (QarraySort(M1int + 1) - QarraySort(M1int))
Else
lowerLimit = QarraySort(1)
End If
If M2 < nofQ - 1 Then
'to be sure that index does not point outside QarraySort
M2int = Int(M2 + 1) 'because the indexing of QarraySort begins from zero
'Interpolation of the upper limit
upperLimit = QarraySort(M2int) + (M2 + 1 - M2int) * (QarraySort(M2int + 1) - QarraySort(M2int))
Else
upperLimit = QarraySort(nofQ)
End If
End Sub 'CalculateConfidenceInterval
Public Function calcB(nYears As Integer, x() As Double, firstYear As Integer, baseYear As Integer, Q As Double) As Double
' calculates the constant B for the equation of linear trend f(t)=Q*t+b.
' The zero point of time axis is the "baseYear"
' Calls the function median
Dim n As Integer 'the number of true values in time series
Dim year As Integer 'the true year of the data value
Dim i As Integer 'index for loop
Dim val() As Double 'array of differences
ReDim val(nYears)
n = 0
For i = 1 To nYears
year = firstYear + i - 1
If x(i) <> MissingValue Then
n = n + 1
val(n) = x(i) - Q * (year - baseYear)
End If
Next i
' the estimate for B is median of the calculated differences
calcB = median(n, val)
End Function ' calcB
Private Function median(nofV As Integer, values() As Double) As Double
' calculates median of values in the array values(), indexed from 1 to nofV
' calls the subroutine sortArray
' is called by the fuction calcB and by the subroutine Sen
Dim i As Integer
Dim sortedValues() As Double
ReDim sortedValues(nofV)
Call SortArray(nofV, values, sortedValues)
If nofV Mod 2 = 0 Then 'nofv is even
i = Int(nofV / 2)
median = (sortedValues(i + 1) + sortedValues(i)) / 2
Else 'nOfvalues is odd
median = sortedValues((nofV + 1) / 2)
End If
End Function 'median
Sub SortArray(ByVal nofV As Integer, values() As Double, sortedValues() As Double)
'This subroutine ranks the values of an array from smallest to largest.
'The sorting method is SELECTION SORT
'The ranked values are stored into a second array called sortedValues.
'Input parameters: nofV - number of values in the array values
' values - values to be ranked, indexed from 1 to nofV
'Subroutine returns the sorted array at sortedValues.
'Is called by the function median and by the subroutine
'CalculateConfidence interval
Dim ind As Integer, i As Integer, j As Integer
Dim minV As Double, maxV As Double
Dim carray() As Double 'the data is first copied to this array
ReDim carray(nofV)
Dim ignoreV As Double 'value that is ignored in carray when sorting
For i = 1 To nofV 'Copy the original array to carray
carray(i) = values(i)
Next i
'Find the smallest and largest value
ind = 1
minV = carray(1) 'initialize the smallest value
maxV = carray(1) 'initialize the largest value
For i = 2 To nofV
If carray(i) < minV Then
minV = carray(i)
ind = i
End If
If carray(i) > maxV Then
maxV = carray(i)
End If
Next i
sortedValues(1) = minV 'the smallest data value
ignoreV = minV - 10 'smaller value than the smallest data value
carray(ind) = ignoreV 'this value is later ignored in sorting
'now sort the values
For j = 2 To nofV
minV = maxV
For i = 1 To nofV
'find the minimum from the rest of the array
If carray(i) <= minV And carray(i) > ignoreV Then
minV = carray(i)
ind = i
End If
Next i
sortedValues(j) = minV
carray(ind) = ignoreV 'from now on this element is ignored
Next j
End Sub 'SortArray
Private Sub fillS()
'Fills the arrays S_nnn of probabilities for two-tailed
' Mann-Kendall test
'The index of tables is the number of data if n=4...10
'Each array entry is an absolute value of the Mann-Kendall
' statistic S, with which the probability that there is no trend
' is less than the probability level p related to the array:
' S_001: p=0.001, S_01: p=0.01, S_05: p=0.05 and S_1: p=0.1.
' Source of values: Gilbert, 1987, Table A18
'Value 9999 indicates that the probability level can not be
' reached with given number of data
Dim n As Integer
For n = 4 To 10
S_001(n) = 9999
S_01(n) = 9999
S_05(n) = 9999
S_1(n) = 9999
Next n
S_001(7) = 21
S_001(8) = 26
S_001(9) = 30
S_001(10) = 35
S_01(6) = 15
S_01(7) = 19
S_01(8) = 22
S_01(9) = 26
S_01(10) = 29
S_05(5) = 10
S_05(6) = 11
S_05(7) = 15
S_05(8) = 18
S_05(9) = 20
S_05(10) = 23
S_1(4) = 6
S_1(5) = 8
S_1(6) = 11
S_1(7) = 13
S_1(8) = 16
S_1(9) = 18
S_1(10) = 21
End Sub
Private Sub Worksheet_SelectionChange(ByVal Target As Range)
End Sub