|
| 1 | +import numpy as np |
| 2 | + |
| 3 | +def OptimalLength(data: np.ndarray): |
| 4 | + """ |
| 5 | + Function calculates the optimal parameter value when using a stationary bootstraping algorithm. |
| 6 | + The method is based on the 2004 paper by Politis & White: |
| 7 | + Dimitris N. Politis & Halbert White (2004) Automatic Block-Length Selection for the Dependent Bootstrap, Econometric Reviews, |
| 8 | + 23:1, 53-70, DOI: 10.1081/ETC-120028836 |
| 9 | +
|
| 10 | + The code was modified compared to Patton's implementation in that it takes as input a one dimensional time-series |
| 11 | + and returns the optimalblock size only for the stationary bootstrap algorithm. |
| 12 | + |
| 13 | + Warning! The minimal size of the time series is 9 elements. |
| 14 | +
|
| 15 | + Example of use: |
| 16 | + >>> import numpy as np |
| 17 | + >>> data = np.array([0.4,0.2,0.1,0.4,0.3,0.1,0.3,0.4,0.2,0.5,0.1,0.2]) |
| 18 | + >>> OptimalLength(data) |
| 19 | + Out[0]: 4.0 |
| 20 | + Args: |
| 21 | + data ... ndarray array containing the time-series that we wish to bootstrap. |
| 22 | + Ex. np.array([-1,0.2,0.3,0.7,0.5,0.1,0.4,0.3,0.5]) |
| 23 | + Returns: |
| 24 | + Bstar ... optimal value of the parameter m Ex. 1 |
| 25 | +
|
| 26 | + Original Matlab version written by: |
| 27 | + James P. LeSage, Dept of Economics |
| 28 | + University of Toledo |
| 29 | + 2801 W. Bancroft St, |
| 30 | + Toledo, OH 43606 |
| 31 | + jpl@jpl.econ.utoledo.edu |
| 32 | +
|
| 33 | + This Python implementation is based on Andrew J. Patton's Matlab code avalible at: |
| 34 | + http://public.econ.duke.edu/~ap172/ |
| 35 | +
|
| 36 | + Implemented by Gregor Fabjan from Qnity Consultants on 12/11/2021. |
| 37 | + """ |
| 38 | + |
| 39 | + n = data.shape[0] |
| 40 | + kn = max(5,np.sqrt(np.log10(n))) |
| 41 | + mmax = int(np.ceil(np.sqrt(n))+kn) |
| 42 | + bmax = np.ceil(min(3*np.sqrt(n),n/3)) |
| 43 | + c = 2 |
| 44 | + |
| 45 | + temp = mlag(data,mmax) |
| 46 | + temp = np.delete(temp,range(mmax),0) |
| 47 | + corcoef= np.zeros(mmax) |
| 48 | + for iCor in range(0,mmax): |
| 49 | + corcoef[iCor] = np.corrcoef(data[mmax:len(data)],temp[:,iCor])[0,1] |
| 50 | + |
| 51 | + temp2 =np.transpose(mlag(corcoef,kn)) |
| 52 | + temp3 = np.zeros((kn,corcoef.shape[0]+1-kn)) |
| 53 | + |
| 54 | + for iRow in range(kn): |
| 55 | + temp3[iRow,:] = np.append(temp2[iRow,kn:corcoef.shape[0]],corcoef[len(corcoef)-kn+iRow-1]) |
| 56 | + |
| 57 | + treshold = abs(temp3) < (c* np.sqrt(np.log10(n)/n)) #Test if coeff bigger than triger |
| 58 | + treshold = np.sum(treshold,axis = 0 ) |
| 59 | + |
| 60 | + count = 0 |
| 61 | + mhat = None |
| 62 | + for x in treshold: |
| 63 | + if (x==kn): |
| 64 | + mhat = count |
| 65 | + break |
| 66 | + count +=1 |
| 67 | + |
| 68 | + if (mhat is None): |
| 69 | + # largest lag that is still significant |
| 70 | + seccrit = corcoef >(c* np.sqrt(np.log10(n)/n)) |
| 71 | + for iLag in range(seccrit.shape[0]-1,0,-1): |
| 72 | + if (seccrit[iLag]): |
| 73 | + mhat = iLag+1 |
| 74 | + break |
| 75 | + if(mhat is None): |
| 76 | + M = 0 |
| 77 | + elif (2*mhat > mmax): |
| 78 | + M = mmax |
| 79 | + else: |
| 80 | + M = 2*mhat |
| 81 | + |
| 82 | + # Computing the inputs to the function for Bstar |
| 83 | + kk = np.arange(-M, M+1, 1) |
| 84 | + |
| 85 | + if (M>0): |
| 86 | + temp = mlag(data,M) |
| 87 | + temp = np.delete(temp,range(M),0) |
| 88 | + temp2 = np.zeros((temp.shape[0],temp.shape[1]+1)) |
| 89 | + for iRow in range(len(data)-M): |
| 90 | + temp2[iRow,:] = np.hstack((data[M+iRow],temp[iRow,:])) |
| 91 | + |
| 92 | + temp2 = np.transpose(temp2) |
| 93 | + temp3 = np.cov(temp2) |
| 94 | + acv = temp3[:,0] |
| 95 | + |
| 96 | + acv2 = np.zeros((len(acv)-1,2)) |
| 97 | + acv2[:,0] = np.transpose(-np.linspace(1,M,M)) |
| 98 | + acv2[:,1] = acv[1:len(acv)] |
| 99 | + |
| 100 | + if acv2.shape[0]>1: |
| 101 | + acv2 =acv2[::-1] |
| 102 | + |
| 103 | + acv3 = np.zeros((acv2.shape[0]+acv.shape[0],1)) |
| 104 | + Counter = 0 |
| 105 | + for iEl in range(acv2.shape[0]): |
| 106 | + acv3[Counter,0] = acv2[iEl,1] |
| 107 | + Counter +=1 |
| 108 | + for iEl in range(acv.shape[0]): |
| 109 | + acv3[Counter,0] = acv[iEl] |
| 110 | + Counter +=1 |
| 111 | + |
| 112 | + Ghat = 0 |
| 113 | + DSBhat = 0 |
| 114 | + LamTemp =lam(kk/M) |
| 115 | + |
| 116 | + for iHat in range(acv3.shape[0]): |
| 117 | + Ghat += LamTemp[iHat]* np.absolute(kk[iHat])*acv3[iHat,0] |
| 118 | + DSBhat += LamTemp[iHat]*acv3[iHat,0] |
| 119 | + DSBhat = 2* np.square(DSBhat) |
| 120 | + |
| 121 | + Bstar = np.power(2*np.square(Ghat)/DSBhat,1/3)*np.power(n,1/3) |
| 122 | + |
| 123 | + if Bstar>bmax: |
| 124 | + Bstar = bmax |
| 125 | + else: |
| 126 | + Bstar = 1 |
| 127 | + return Bstar |
| 128 | + |
| 129 | + |
| 130 | +def mlag(x: np.ndarray,n)-> np.ndarray: |
| 131 | + """ |
| 132 | + Returns a numpy array in which the k-th column is the series x pushed down (lagged) by k places. |
| 133 | + |
| 134 | + Example of use |
| 135 | + >>> import numpy as np |
| 136 | + >>> x = np.array([1,2,3,4]) |
| 137 | + >>> n = 2 |
| 138 | + >>> mlag(x,n) |
| 139 | + Out[0]: array([[0, 0], |
| 140 | + [1, 0], |
| 141 | + [2, 1], |
| 142 | + [3, 2]]) |
| 143 | + The function was tested passing a numpy array (ndarray) as input and requires numpy to run. |
| 144 | + Args: |
| 145 | + x ... ndarray array for which the lagged matrix is calculated. np.array([1,2,3,4]) |
| 146 | + n ... integer specifying how many lags does the function consider |
| 147 | + Returns: |
| 148 | + xlag... ndarray contining the k-th lagged values in the k-th column of the matrix |
| 149 | +
|
| 150 | + Original Matlab version written by: |
| 151 | + James P. LeSage, Dept of Economics |
| 152 | + University of Toledo |
| 153 | + 2801 W. Bancroft St, |
| 154 | + Toledo, OH 43606 |
| 155 | + jpl@jpl.econ.utoledo.edu |
| 156 | +
|
| 157 | + This Python implementation is based on Andrew J. Patton's Matlab code avalible at: |
| 158 | + http://public.econ.duke.edu/~ap172/ |
| 159 | + |
| 160 | + Implemented by Gregor Fabjan from Qnity Consultants on 12/11/2021 |
| 161 | + """ |
| 162 | + nobs = x.shape[0] |
| 163 | + out = np.zeros((nobs,n)) |
| 164 | + for iLag in range(1,n+1): |
| 165 | + for iCol in range(nobs-iLag): |
| 166 | + out[iCol+iLag,iLag-1] = x[iCol]; |
| 167 | + return out |
| 168 | + |
| 169 | + |
| 170 | +def lam(x: np.ndarray)-> np.ndarray: |
| 171 | + """ |
| 172 | + Returns the value at points x of the Trapezoidal function. Trapezoidal funcion maps all numbers bigger than 1 or smaller than -1 to zero. |
| 173 | + Values between -1/2 to 1/2 to 1 and the rest either on the line connecting (-1,0) to (-1/2,1) or (1/2,1) to (1,0). |
| 174 | +
|
| 175 | + Example of use: |
| 176 | + >>> import numpy as np |
| 177 | + >>> x = np.array([0.55]) |
| 178 | + >>> lam(x) |
| 179 | + Out[0]: array([0.9]) |
| 180 | +
|
| 181 | + Args: |
| 182 | + x ... ndarray array of points on which we wish to apply the trapezoidal mapping. |
| 183 | + Ex. np.array([-1,0.2,0.3,0.7]) |
| 184 | + Returns: |
| 185 | + out ... ndarray of mapped points Ex. array([0. , 1. , 1. , 0.6]) |
| 186 | +
|
| 187 | + Original Matlab version written by: |
| 188 | + James P. LeSage, Dept of Economics |
| 189 | + University of Toledo |
| 190 | + 2801 W. Bancroft St, |
| 191 | + Toledo, OH 43606 |
| 192 | + jpl@jpl.econ.utoledo.edu |
| 193 | +
|
| 194 | + This Python implementation is based on Andrew J. Patton's Matlab code avalible at: |
| 195 | + http://public.econ.duke.edu/~ap172/ |
| 196 | + |
| 197 | + Implemented by Gregor Fabjan from Qnity Consultants on 12/11/2021. |
| 198 | + """ |
| 199 | + |
| 200 | + nrow = x.shape[0] |
| 201 | + out = np.zeros(nrow) |
| 202 | + for row in range(nrow): |
| 203 | + out[row] = (abs(x[row])>=0) * (abs(x[row])<0.5) + 2 * (1-abs(x[row])) * (abs(x[row])>=0.5) * (abs(x[row])<=1) |
| 204 | + return out |
| 205 | + |
| 206 | + |
0 commit comments