Quantitative processes are necessarily complex. Explaining them requires a deep understanding of the beliefs. Confirming those beliefs requires a deep understanding of the maths. It is a snowy, blistering cold day in London. A perfect time to be introspective. So with your consent, I’d like to review a few methods for confirming factor premia – tests for monotonicity.

There are sound theoretical reasons for expecting cross sectional differences in exposures to result in monotonic difference in returns. If one constructs n-tile portfolios, the lowest n-tile should under-perform the second n-tile and so on with the highest n-tile providing the highest performance.

The tests to confirm this outcome are necessarily complex. The first of these test is by Andrew Patton and Allen Timmerman. Another by Joseph P Romano and Michael Wolf. Both examine the differences in adjacent n-tile portfolios. How often do the returns of portfolio n-tile_{k-1} outperform n-tile_{k}? That should not happen too often. Patton and Timmerman asks if the relationship is monotonic increasing. Romano and Wolf respond with not so fast. There are cases where some n-tiles are larger than the previous but also larger than the next n-tile. Not strictly increasing. When this is the case, are the violations small enough in occurrence to still accept monotonicity?

The difference is that the Patton and Timmerman test allows for a wider class of outcomes to pass the monotonicity test. Romano Wolf is stricter. From the conclusion section of Romano and Wolf:

In a recent paper, Patton and Timmermann (2010) propose a new test, taking all categories of the underlying characteristic into account. Compared to previous related proposals, they are, to our knowledge, the first ones to postulate a strictly increasing relation as the alternative hypothesis of the test, rather than as the null hypothesis. This is the correct formulation if the goal is to establish strict monotonicity with a quantifiable statistical significance. In addition, they postulate a weakly decreasing relation as the null hypothesis. Compared to allowing, more generally, also for a non-monotonic or weakly increasing relation under the null, this approach results in a smaller critical value and thereby in higher power of the test. But if a non-monotonic or weakly increasing relation is not ruled out, the test can break down and falsely ‘establish’ a strictly increasing relation with high probability. Therefore, the test of Patton and Timmermann (2010) should only be used as a test for the direction of (existing) monotonicity.

In this paper, we have proposed some alternative tests that also allow for a non-monotonic or weakly increasing relation under the null and that successfully control the probability of a type 1 error. These tests do not require modeling the dependence structure of the time series nor the covariance matrix of the observed return differentials. A suitable bootstrap method is used to implicitly account for these unknown features of the underlying data generating mechanism. As is unavoidable, such tests have lower power compared to the test of Patton and Timmermann (2010). This is the price one has to pay to safe-guard against the possibility of falsely ‘establishing’ strict monotonicity (beyond the nominal significance level of the test) in general settings.

Here is my code. I’ve altered Patton and Timmerman’s block bootstrap code with a version inspired by Kevin Sheppard’s code, it is much faster. Patton’s code is quite clear to understand, but not fast.

Change to :

% generating the time indices for the bootstrapped data:

bootdates = block_bootstrap(T,bootreps,block_length,rand_state);

%bootdates = stat_bootstrap_function_21(T,bootreps,block_length,rand_state); old code

— the disclaimer bits and BSD 2-clause licence—-

** If you injure yourself using any of this, you’ve only yourself to blame.**

** Copyright 2018 louis scott lscott@kiemaadvisors.com**

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

----- the demo: demo_mono_rw.m ---- % --------------------------------------------------------- % AUTHOR: lscott@kiemaadvisors.com % --------------------------------------------------------- % % Copyright 2018 louis scott % This file is released under the BSD 2-clause license. % The BSD 2-clause license is located in the license directory, % or type licence at the matlab/octave command prompt % % --------------------------------------------------------- % for a patton and timmerman test, substitute % out_pt = MR_test_22(data,pm.bootreps,0,pm.bootreps,pm.rand_state) pm.kernel = 'quad'; pm.bootreps = 2000; pm.rand_state = 123456; data = randn(2000,5); for ii = 1:3 rescale = .0225 +(ii-1)*.004; datam = data + repmat(rescale*[1:5],2000,1); out_rw = MR_test_RW(datam,pm); disp(out_rw(4)); eval(['subplot(31' num2str(ii) ');']); semilogy(cumprod(1+.01*datam),'Linewidth',2); axis([1 2000 .8 100]) legend({'1','2','3','4','5'},'Location','NorthWest') title(sprintf('Probability of non-monotonicity: %4.3f',out_rw(4))) end print('demo_mono_rw','-dpng') ----- Romano wolf in matlab/octave ----- function res = MR_test_RW(data,params) % % Function to compute the "monotonic relationship" tests in: % Joseph P Romano and Michael Wolf, Testing for Monotonicity in Expected Asset Returns % https://papers.ssrn.com/sol3/papers.cfm?abstract_id=1858761 % % --------------------------------------------------------- % AUTHOR: lscott@kiemaadvisors.com % --------------------------------------------------------- % % Copyright 2018 louis scott % This file is released under the BSD 2-clause license. % The BSD 2-clause license is located in the license directory, % or type licence at the matlab/octave command prompt % % --------------------------------------------------------- [T,n] = size(data); res = nan(4,1); params = setdefaults(params); if params.direction==0; params.direction = sign(mean(data(:,n)-data(:,1))); end if abs(params.direction)<2 % then need to difference the data diffdata = data(:,2:end)-data(:,1:end-1); % differences across the columns of the data diffdata = params.direction*diffdata; % changing the sign of these if want to look for a *decreasing* pattern rather than an *increasing* pattern else diffdata = data; % data are already in differences diffdata = sgn(params.direction)*diffdata; % changing the sign of these if want to look for a *decreasing* pattern rather than an *increasing* pattern n = n + 1; % if data are already differenced, then this is like having one more column in the raw data end dmuhat = mean(diffdata)'; if abs(params.direction)<2 % if data are already in differences, then I cannot do the t-test. this must be done on the actual returns % now getting the t-statistic and p-value from the usual t-test temp = nwest(params.direction*(data(:,end)-data(:,1)), ones(T,1), floor(4*((T/100)^(2/9)))); % truncation lag for Newey-West standard errors, see Newey and West (1987, 1994) res(1) = temp.tstat; res(2) = normcdf(-temp.tstat); end % generating the time indices for the bootstrapped data: r_{i}^{*} bootdates = block_bootstrap(T,params); temp = nan(params.bootreps,n-1); tempS = temp; tstar = temp; tempSE = temp; sigma_hatstar = nan(n-1,1); % Algorithm 4.1 % Generate bootstrap data {r_{1}^{*}, r_{2}^{*} , . . . , r_{T}^{*} } and % compute statistics \hat{∆}_{T,i} and \hat{σ}_{T,i} % from the data, for i = {1,\dots,N} % patton and timmerman results for later use in creating conservative null % \delta_{0,i} res_PT = MR_test_22(data,params.bootreps,params.direction,params.block_length,params.rand_state); del_cons = res_PT(5,2); % compute \hat{\delta} and t_{i}^{*} for i = {1,\dots.N} % % HAC standard errors are computed with % Jochen Heberle and Cristina Sattarhoff, % A Fast Algorithm for the Computation of HAC Covariance Matrix Estimators for ii=1:n-1 temp2 = diffdata(:,ii); delta_hat = mean(temp2(bootdates))'; sigma_hatstar(ii) = sqrt((1/T)*HAC_hs(delta_hat, params.kernel, params.block_length)); for jj=1:params.bootreps tempSE(jj,ii) = sqrt((1/T)*HAC_hs(temp2(bootdates(:,jj)), params.kernel, params.block_length)); end %sigma_hatstar = mean(tempSE(bootdates))'; % Compute t_{T}^{min,*} \equiv min_{i}t_{T,i}^{*} % Repeat this process B times, resulting in statistics % t_{T,1}^{min,*}, t_{T,2}^{min,*} \dots t_{T,B}^{min,*} temp(:,ii) = delta_hat - dmuhat(ii) + del_cons; % the mean of each of the bootstrap shuffles of the original data, % minus the mean of the original data (the re-centering bit) tempS(:,ii) = (delta_hat - dmuhat(ii) + del_cons)./tempSE(ii); % t_{i}^{*} for i = {1,\dots.N} eq 4.1 end tempm = min(temp,[],2); % looking at the minimum difference in portfolio means, for each of the bootstrapped data sets Jstat = min(dmuhat); % the test statistic res(3,1) = mean(tempm>Jstat); % the p-value associated with the test statistic % t_{T,1}^{min,*}, t_{T,2}^{min,*} \dots t_{T,B}^{min,*} tempSm = min(tempS,[],2); % looking at the minimum STUDENTISED difference in portfolio means, for each of the bootstrapped data sets JstatS = min(dmuhat./sigma_hatstar); % the STUDENTISED test statistic res(3) = JstatS; res(4) = mean(tempSm>JstatS); waithere = true; % Based on remark 4.2 alternative way to implement the Cons test. end function params = setdefaults(params); if(~isfield(params,'bootreps')) params.bootreps = 1000; end if(~isfield(params,'direction')) params.direction = 0; end if(~isfield(params,'block_length')) params.block_length = 10; end if(~isfield(params,'kernel')) params.kernel = 'Parzen'; end if(~isfield(params,'rand_state')) params.rand_state = []; end if(isempty(params.bootreps)) params.bootreps = 1000; end if(isempty(params.direction)) params.direction = 0; end if(isempty(params.block_length)) params.block_length = 10; end if(isempty(params.kernel)) params.kernel = 'Parzen'; end if(isempty(params.rand_state)) rand('state',sum(100*clock)); params.rand_state = rand("state"); end if(~isnumeric(params.bootreps)) error('bootreps must be a number') end if(~isnumeric(params.block_length)) error('block_length must be a number') end if(~isnumeric(params.direction)) error('direction must be a number') end if(params.block_length< 0) error('block_length must be > 0') end if(params.bootreps < 0) error('bootreps must be > 0') end if(params.block_length < 0) error('block_length must be > 0') end end ----- Patton and Timmerman available here: ----- http://public.econ.duke.edu/~ap172/Patton_Timmermann_MR_test_toolbox.zip ----- The Block bootstrap ----- % % --------------------------------------------------------- % AUTHOR: lscott@kiemaadvisors.com % --------------------------------------------------------- % % Copyright 2018 louis scott % This file is released under the BSD 2-clause license. % The BSD 2-clause license is located in the license directory, % or type licence at the matlab/octave command prompt % % --------------------------------------------------------- % Notes: A variant on the code by Kevin Sheppard kevin.sheppard@economics.ox.ac.uk % see https://www.kevinsheppard.com/UCSD_GARCH % his is released under the BSD licence, so this is compliant. % function indices = block_bootstrap(t,br,w,rand_state) if(isstruct(br)) w = br.block_length; rand_state = br.rand_state; tt= br.bootreps; br = tt; end rand('state',rand_state); % Compute the number of blocks needed s=ceil(t/w); % Generate the starting points Bs=floor(rand(s,br)*t)+1; indices=zeros(s*w,br); index=1; % Adder is a variable that needs to be added each loop adder=repmat((0:w-1)',1,br); for i=1:w:t indices(i:(i+w-1),:)=repmat(Bs(index,:),w,1)+adder; index=index+1; end ndx = (indices > t); indices(ndx) = indices(ndx) - t; end ----- fast covariance estimation ------ see Jochen Heberle and Cristina Sattarhoff, A Fast Algorithm for the Computation of HAC Covariance Matrix Estimators, Here is my implementation: function res = HAC_hs(A, kernel,bw) % % try to initialize fftw for speed gains % see eqns 38-39 % Jochen Heberle and Cristina Sattarhoff, A Fast Algorithm for the Computation of HAC Covariance Matrix Estimators, % Econometrics 2017, 5, 9, pp % % --------------------------------------------------------- % AUTHOR: lscott@kiemaadvisors.com % --------------------------------------------------------- % % Copyright 2018 louis scott % This file is released under the BSD 2-clause license. % The BSD 2-clause license is located in the license directory, % or type licence at the matlab/octave command prompt % % --------------------------------------------------------- % Notes: my translation of their R code into matlab/octave % [N,q] = size(A); wstar = HACwts(kernel,N,bw); % Compute the eigenvalues $ λ i (i = 1, . . . , 2N) of C(ω^∗) $ using Equation (34) with wstar = [1; wstar(1:N-1); 0 ; wstar(N-1:-1:1)]; wstar = real(fft(wstar)); % Construct the matrix $A^∗$ with dimension 2N × q from FF = [A; zeros(N,q)]; % Determine $V^∗ A^{∗}_{j} $ given by the DFT of $A^{∗}_{j}$. FF = fft(FF); % Multiply for all $i ∈ { 1, . . . , 2N }$ the i-th entry of the vector $V^{*} A^{∗}_{j}$ % with the eigenvalue $λ_{i}$ , in order to construct $ΛV^{*}A^{∗}_{j}$ . FF = FF.*repmat(wstar,1,q); % Determine $C(ω^∗) A^{∗}_{j} = VΛV ∗ A^{∗}_{j} given by the IDFT of ΛV^{*}A^{∗}_{j}$ . FF = real(ifft(FF)/(2*N)); FF = FF(1:N,:); % $\hat{S}_N = \frac{1}{N}{A}'T(ω)A. where FF = T(ω)A.$ res = (1/N)*(A'*FF); end ------HAC_hs relies on the code HAC_wts ----- function ww = HACwts( kernel, dN, bw) % % with kernel = {"Truncated", "Bartlett", "Parzen", "Tukey-Hanning", "Quadratic Spectral"} % bw the bandwidth of nonzero weights % dN matrix dimension that calls HAC % % --------------------------------------------------------- % AUTHOR: lscott@kiemaadvisors.com % --------------------------------------------------------- % % Copyright 2018 louis scott % This file is released under the BSD 2-clause license. % The BSD 2-clause license is located in the license directory, % or type licence at the matlab/octave command prompt % % --------------------------------------------------------- % Notes: my translation of their R code into matlab/octave % Jochen Heberle and Cristina Sattarhoff, A Fast Algorithm for the Computation of HAC Covariance Matrix Estimators, % Econometrics 2017, 5, 9, pp % ww = zeros(dN,1); switch lower(kernel(1:4)) case "trun" ww(1:bw) = 1; case "bart" c1 = (bw - 1)/2; ww(1:bw) = 1 - ((1:bw)/(bw + 1)); case "parz" c1 = (1:floor(bw/2))/bw; c2 = ((floor(bw/2)+1):bw)/bw; ww(1:length(c1)) = 1 - 6*c1.^2 + 6*c1.^3; ww(length(c1)+1:bw) = 2*(1-c2).^3; case "tuke" ww(1:bw) = (1 + cos(pi*((1:bw)/bw)))/2; case "quad" aa = pi*((1:dN)'/bw)/5.0; ww = 1.0./(12*aa.^2).*(sin(6*aa )./(6*aa) - cos(6*aa)); otherwise error('kernel case not understood'); end end --- licence.m function licence() display('Copyright 2018 louis scott lscott@kiemaadvisors.com') disp('This file is released under the BSD 2-clause license.') disp('') disp('Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:') disp('') disp('1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.') disp('') disp('2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.') disp('') disp('THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.') end