## Factor diagnostics: monotonicity with Romano and Wolf, Patton and Timmerman 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.

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  ----

% ---------------------------------------------------------

% ---------------------------------------------------------
%
%  This file is released under the BSD 2-clause license.
%  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.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
%

% ---------------------------------------------------------

% ---------------------------------------------------------
%
%  This file is released under the BSD 2-clause license.
%  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 -----
%
% ---------------------------------------------------------

% ---------------------------------------------------------
%
%  This file is released under the BSD 2-clause license.
%  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
for i=1:w:t
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
%
% ---------------------------------------------------------

% ---------------------------------------------------------
%
%  This file is released under the BSD 2-clause license.
%  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
%
% ---------------------------------------------------------

% ---------------------------------------------------------
%
%  This file is released under the BSD 2-clause license.
%  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;
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()

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



## The NAG libraries in Julia.

NAG.jl  The NAG libraries in Julia

Friends, it is time to get on board the Julia train.  With the release of v0.6 and the promise of v1.0 around the corner, it seemed an opportune time to dig in.  Hating manuals, I decided to learn by contributing to a package in need of a developer, the NAG library.   The integration is discussed in a post at the NAG blog by NAG’s Reid Atcheson and Julia’s Andy Greenwell.  Reid gave a recent workshop on the use of NAG on the Intel Phi,  Andy gave a great talk at JuliaCon on C,C++ integration.   Both well worth your time.  Julia’s promise of being the one language for high performance computing and scientific environment is what attracted me, and the use of the NAG libraries certainly is a useful confluence of these two areas.

So learning Julia. A joy really.  Have a look at the JuliaCon 2016 workshops on youtube.  In particular, David Saunders Intermediate Julia workshop, Arch Robison’s high performance Julia workshop, Tom Sargent’s keynote on the use of Julia in macro economics , and the FRB’s Erika Moszkowski’s talk on a specific econometric topic DSGE models in Julia.  I could go on… I watched hours of the conference, the level of quality here is stellar.

With that background, I began coding by getting to the existing code NAG.jl on github, cloning it and beginning.  I also took a look at GSL.jl, the wrapper for the GNU scientific library coded by Yichao Yu to get a sense of good practice for something this large.  What I found is a reliance on the metaprogramming style of writing Julia.  Julia, like lisp is homoiconic, the code looks like its abstract syntax tree.  That makes it easier to parse and extend.  There are two packages that take this on, clang.jl giving access to a large subset of llvm clang, and BaseTestAuto.jl for automatted testing.   Now step two could not be learn clang.jl – its too big a fish to reel in as I learn.  Instead I recognised that a number of special functions in chapter S take real or complex numbers and return real or complex numbers.  These simple file signatures in the nags.h header file did not require clang to automate the wrapping of these functions.  I could rely on a little bash with sed to grab what was needed.  Here is what that looks like:

fn_type_str = "return complex, inputs (complex z, NagError *fail)"
s_dict = Dict{Symbol, Symbol}(:s14agc => :nag_complex_log_gamma,:s15ddc => :nag_complex_erfc)
for nagf in (:s14agc, :s15ddc,
)
@eval begin
@doc "$(section_str) : Fns that$(fn_type_str) $(nag_doc_str) "$nagf
export $(s_dict[nagf]) function$(s_dict[nagf])(z::Number)
reset_nag_error()

## Yesterday, risk was such an easy game to play.

There is no need for nuance in explaining yesterdays market – we have an appetite for risk. The move in equities was as broad as one is ever likely to see, with all sectors in line with their ex-ante beta. European equities being the laggard. Bonds also followed suit, with the high yield corporate ETF HYG (+.9%) and junk ETF JNK (+.8%) having a strong day and treasuries down, the 20 year TLT ETF falling (-1.3%).

Commodities fared about as well as corporates, with the notable exception of Natural gas, UNG fell sharply on the day (-3.5%). Note that UNG has a high beta with the AGG, Carry and TIPs, and it was not a good day to be associated with TIPs.

By far the simplest story here is this: The outliers being Natural gas, and near and mid term VIX, the graph clusters by asset class above the outliers. It was not a good day to be sensitive to the AGG. Precious, base metals, natural gas and oil are sensitive to inflation beta – there was quite the spead in performance among these.  The Natural gas coefficient is not significant despite being quite large. Oil, emerging markets, and Asia ex Japan have large and significant loadings on the returns to the carry trade.  These are followed by commodities in general.

## Beating Benchmarks with No Skill*

##### As in Lochoff 1998, the benchmark is the Barclay (Lehman) Aggregate bond index.   Where Lochoff used the historical means and returns to construct a duration matched portfolio that shorts the risk free and goes long (with leverage) on short duration assets – say from three months to 3 years.   I also use historical data but in contrast employ an S-shaped utility function to state my relative preference for a return distribution.  We use Utility =I(x>0)* x^m- L*I(x<0)*(-x)^m with (L,m) = (2,.95), (2,.8) and (3,.95) and feed it to our FSO.  That is we have an asymmetric preference, we dislike a negative outcome 2-3 times the magnitude of a positive outcome.   Distributions with negative tails then need considerable mass on the positive side to compensate for the tails.  See Cremers, Kritzman and Page 2004, Adler and Kritzman 2005 and Hagströmer 2007 for details of FSO and the S-shaped utility function.  We then equal weight each of the 3 allocations.
                           Downside  Return/    Return/
Mean      Stdev    Stdev     Stdev     Downside Stdev
AGG     0.0600    0.0657    0.0580    0.9131    1.0337
Port    0.0829    0.0747    0.0513    1.1089    1.6167

Portfolio Beta: .871

## Risk Appetite Indicator

##### The Indicator is makes use of the construction methodology in Stock and Watson 2004,  Combination forecasts of output growth in a seven country data set, as well as the more recent paper by Rapach, Strauss, and Zhou, Out-of-sample equity premium prediction: Combination forecasts and links to the real economy and a working paper Ibisevic, Snijder, Lai, and Zhang, Combining forecasts: Improved approach and analysis of various frameworks

The distribution of QRAI exhibits a tri-modality.
Peaks corresponding to a normal regime, a quiet regime and a bull run regime.

Risk Appetite Index Results

The Index is used to enter and exit the S&P500 when the level of the index is low. The Anatolyev & Gerko test of Expected Profitability is employed as well as a related regression that shows the profitability of both the long and short legs of the strategy.

Anatolyev Gerko EP Statistic

Ordinary Least-squares Estimates
Dependent Variable =           SP500
R-squared      =    0.0106
Rbar-squared   =    0.0101
sigma^2        =    0.0002
Durbin-Watson  =    2.1638
Nobs, Nvars    =   2103,     2
***************************************************************
Variable      Coefficient      t-statistic    t-probability
const            0.002848         4.599309         0.000005
qri EP           0.002938         4.743630         0.000002

## Full Scale Optimization with REITs

##### The distributions of the three assets show the large tails in the REITs and the relative safety of the Ten year. We use the S-shaped utility function as our canonical utility function. It is asymmetric and defined as follows. For returns less than a threshold – zero in our case – the utility is 1.5*Loss^.95 for gains it is Gains^.95. Here is an image of the S-shaped utility for each allocation of (REITs, 10 year, SP500) when the exponent is .5. Quite difficult to do gradient descent. The S-shaped Utility for our problem when exponent is .5 presents an issue to traditional gradient descent.

##### The MV optimization use a RAP of .5 on annualized returns. We allow the covariance to be based on an expanding window, and an exponentially weighted moving average with lambda = 99.5. The portfolio construction limits the maximum weight to 85% of the portfolio, and the weights sum to one. The time series of allocations for all four portfolios differs dramatically across MV and FSO

##### The optimization RAP and the S-shaped asymmetry were chosen so that the two approaches yielded similar realized portfolio returns. Here we show the portfolio performance for our four cases, FSO and MV with expanding (F for full) and ewma weighted (V for varying) windows. The annualized returns are comparable, the downside variance is slightly better for the FSO’s and consequently the Sortinos are better. We also show the 1/N portfolio. Portfolio returns for each way of constructing portfolios gives the advantage to FSO

##### While the performance differences are small the tails are different. Here I show the realized average shortfall below -1%. The FSO distributions have smaller left tails and this is significant at the 2% level for the Kolomogorov- Smirnov distance test. The fatter tails in the MV optimizations are significantly different from their FSO counterpart using a KS test. Never mind the 1/N which isn't in the same ballpark.

## Toeplitz matrices as families of shrinkage covariance estimators

##### My current project requires that I examine the risk of a large number of portfolios so I need to create some sensible estimates of risk.  Which in part, means looking for robust shrinkage estimators of covariance matrices.  A requirement is that the resulting matrices are stable to random perturbations.  Specifically if I add a little random noise, I should still get nearly the same answer. So I turn to Toeplitz matrices, with some cool results.To see which matrices and transformations of matrices are stable, it helps to start with some very structured matrices that look like risk covariances, and therefore allow us to examine if all that structure is preserved when we kick them with some random data. We use Toeplitz matrices – like risk model covariance matrices, they are symmetric can be shown to have real eigenvalues and in some special cases, their eigenvectors have closed form solutions.It is also known that certain Toeplitz matrices can be proven to have eigenvectors that are periodic.So I start with a very simple diagonal Toeplitz matrixNth=2^9; bx=[2*Nth:-2:1]/(2*Nth); vcv=toeplitz(bx.^4); image(vcv)The eigenvectors have a lot of structure, this is the first image in this post above. The eigenvalues are unremarkable and not the story here. For the record, the eigenvalues are real and well behaved in our examples.[ev,ed] = eig(vcv); image(ev); Even if we preturb the matrix with random data vcvDr= vcv + 1e-8*rand(Nth,Nth); [evDr,edDr] = eig(vcvDr); vcvDr2= vcv + 1e-3*rand(Nth,Nth); [evDr2,edDr2] = eig(vcvDr2); Symmetric Toeplitz matrix The first five eigenvectors are well behaved under perturbationThe first of four plot shows the first 5 eigenvectors of the original matrix, then to the right are the eigenvectors of a disturbed matrix vcvDr where we add a random perturbation on the order of 10^-8. The bottom left is another larger random perturbation of size 10^-3. Finally we show the difference in absolute value for the entries in the first five eigenvectors. We use absolute value as the eigenvectors are unique up to sign differences. No problem with stability here. What is notable is that this is similar to Ledoit-Wolf where the diagonal is also constant – they use a market weighted average variance.  The off-diagonals in Ledoit-Wolf are average covariance. Here the off-diagonals are not constant but allow for the linkages between assets to fall off as matrix polynomial – here we just show the 4th power. So the matrix allows for some assets to be tightly coupled and others to be decoupled.  In fact the design lends itself to be decomposed into blocks where the blocks could be asset class covariances. Those of you who have dealt with fixed income risk models know that this tightly coupled group is difficult to work with. One answer is to only use a few say three principal components.  All good and fine, what I am addressing allows for that but more to the point, the stability issue can hit the largest eignevalues and vectors. Another reason to consider the stability question is that the use of a toeplitz matrix with noise allows you to explore the structure that is preserved as you add random noise. When we allow for the Toeplitz design to be relaxed and allow for the main diagonal to vary exponentially – that is allowing for a parametrized variance where some assets are substantially more volatile than others, then we need to ask about stability. First with an exponent of .995:  for ith= 2:Nth for jth= ith:Nth vcvD(ith,jth) = xpwt*vcvD(ith-1,jth-1); vcvD(jth,ith) = vcvD(ith,jth); end end figure; imagesc(vcvD)The eigenvalues are structured, but now their is an exponential decay together with the periodicity.  We also see that the addition of noise leaves the eigenvalues stable. The symmetric weighted exponential matrix (.995) Structure of the eigenvectors for the symmetric weighted exponential matrix (.995) The eigenvalue structure for the first five ev's. (lambda=.995)When we allow for an exponent of .985, the decay is large and the matrix looks quite different.The eigenvalues however are still quite similar.  So are the perturbed eigenvectors, but the noise is apparent and of the same scale as the vector coefficients. If you are thinking of using shrinkage estimates that vary along the diagonal you should be mindful of this result.

The matlab code for this:
Nth=2^9; bx=[2*Nth:-2:1]/(2*Nth); for xpwt = [1.0 .995 .985] vcv=toeplitz(bx.^4); [ev,ed] = eig(vcv); vcvD = vcv; for ith= 2:Nth for jth= ith:Nth vcvD(ith,jth) = xpwt*vcvD(ith-1,jth-1); vcvD(jth,ith) = vcvD(ith,jth); end end figure; imagesc(vcvD)

 % [evD,edD] = eig(vcvD); vcvDr = vcvD + 1e-8*rand(Nth,Nth); [evDr,edDr] = eig(vcvDr); vcvDr2= vcvD + 1e-3*rand(Nth,Nth); [evDr2,edDr2] = eig(vcvDr2); [ntoss,ndx]=sort(diag(edD),'descend'); [ntoss,ndxr]=sort(diag(edDr),'descend'); [ntoss,ndxr2]=sort(diag(edDr2),'descend'); figure; imagesc(min(max(evD(:,ndx),-.07),.07)) 

figure subplot(221); plot(evD(:,ndx(1:5))); axis([1 Nth -.07 .07]) subplot(222); plot(evDr(:,ndxr(1:5))); axis([1 Nth -.07 .07]) subplot(223); plot(evDr2(:,ndxr2(1:5))); axis([1 Nth -.07 .07]) subplot(224); plot([abs(evD(:,ndx(1:5))) - abs(evDr2(:,ndxr2(1:5)))]); axis([1 Nth -.07 .07]) end

##### References A. Bottcher and S.M. Grudsky, Toeplitz Matrices, Asymptotic Linear Algebra, and Functional Analysis, Birkhauser, 2000. A. Bottcher and B. Silbermann, Introduction to Large Truncated Toeplitz Matrices, Springer, New York, 1999. W. Cheney, Introduction to Approximation theory, McGraw-Hill, 1966. T. A. Cover and J. A. Thomas, Elements of Information Theory, Wiley, New York, 1991. P. J. Davis, Circulant Matrices, Wiley-Interscience, NY, 1979. D. Fasino and P. Tilli, “Spectral clustering properties of block multilevel Hankel matrices, Linear Algebra and its Applications, Vol. 306, pp. 155–163, 2000. F.R. Gantmacher, The Theory of Matrices, Chelsea Publishing Co., NY 1960. R.M. Gray Toeplitz and Circulant Matrices: A review. 206, Now Publishers R. M. Gray, “On the asymptotic eigenvalue distribution of Toeplitz matrices,” IEEE Transactions on Information Theory, Vol. 18, November 1972, pp. 725–730. Ulf Grenander and Gabor Szego, Toeplitz forms and their application, American Mathematical Soc., 1984 Lloyd Nick Trefethen and David Bau, Numerical Linear Algebra, SIAM 1997 W. F. Trench, “Asymptotic distribution of the even and odd spectra of real symmetric Toeplitz matrices,” Linear Algebra Appl., Vol. 302-303, pp. 155–162, 1999. W. F. Trench, “Absolute equal distribution of the spectra of Hermitian matrices,” Lin. Alg. Appl., 366 (2003), 417–431.
Posted in Code, matlab, Matrix Algebra | Tagged , , | 1 Comment