Price discovery in high resolution: 
computational appendix

Joel Hasbrouck*
October 1, 2018

*Department of Finance
Stern School of Business, NYU
44 West 4th St.
New York, NY 10012
jhasbrou@stern.nyu.edu

Abstract

Estimating and forecasting long-lag high-resolution multivariate time series models pose computational challenges. This appendix describes several techniques used in the paper. Most are based on standard polynomial distributed lag (PDL) and sparse matrix methods.  I discuss these in simple univariate and bivariate settings, along with extensions relevant to the problem at hand.

This document includes Mathematica code at the end. The Matlab software uses additional conventions and will be described in a separate document. There are two versions of this appendix: a pdf document (with a 'pdf' extension) and a Mathematica notebook (with an 'nb' extension). The latter is executable in Mathematica.

A review of polynomial distributed lags (PDLs)

For bivariate time series {(Subscript[x, t]),(Subscript[y, t])}, consider a standard linear lagged regression model:

(Subscript[y, t])=(Subscript[β, 1])(Subscript[x, t-1])+(Subscript[β, 2])(Subscript[x, t-2])+...+(Subscript[β, ℒ])(Subscript[x, t-ℒ])+(Subscript[e, t])

If (Subscript[x, t])=(Subscript[y, t]) this is an autoregression. The data sample is {(Subscript[x, t]),(Subscript[y, t]):t=1,...,𝒯}. The values can be arranged for estimation as y=X β+e, where y=(({{(Subscript[y, 1]), (Subscript[y, 2]),..., (Subscript[y, 𝒯]))}})ᵀ, e is a 𝒯×1 vector of errors, and X is the lagged data matrix, that is a 𝒯×ℒ matrix of the lagged xs:

X=(({{(Subscript[x, 0]), (Subscript[x, -1]), (Subscript[x, -2]), ⋯, (Subscript[x, -ℒ+1])}, {(Subscript[x, 1]), (Subscript[x, 0]), (Subscript[x, -1]), …, (Subscript[x, -ℒ])}, { , ⋮,  ,  ,  }, {(Subscript[x, t]), (Subscript[x, t-1]), (Subscript[x, t-2]), ⋯, (Subscript[x, t-ℒ+1])}, { , ⋮,  ,  ,  }, {(Subscript[x, 𝒯-1]), (Subscript[x, 𝒯-1]), (Subscript[x, 𝒯-1]), ⋯, (Subscript[x, 𝒯-ℒ+1])}}))

Here, values of (Subscript[x, t]) for t<1 can be handled in various ways. If E x=0 we can set the pre-sample values to zero. Alternatively, we can discard rows that contain any pre-sample values; y is adjusted to be conformable. The OLS estimates of β are (OverHat[β])=(((XᵀX))^(-1))Xᵀy.
In the autoregressive case, y may be forecast subsequent to an initial disturbance (Subscript[e, 0]) by recursion. The forecasts, denoted (Subscript[y, t]^( * )), are:

(Subscript[y, 1]^( * ))=(Subscript[e, 0])
(Subscript[y, 2]^( * ))=(Subscript[β, 1])(Subscript[y, 1]^( * ))
(Subscript[y, 3 ]^( * ))=(Subscript[β, 1])(Subscript[y, 2]^( * ))+(Subscript[β, 2])(Subscript[y, 1]^( * ))
⋮
(Subscript[y, t]^( * ))=(Subscript[β, 1])(Subscript[y, t-1]^( * ))+(Subscript[β, 2])(Subscript[y, t-2]^( * ))+...+(Subscript[β, ℒ])(Subscript[y, t-ℒ]^( * ))

This is the same as the lagged regression model, but since we are forecasting, we are recursively calculating the system assuming that all disturbances, except (Subscript[e, 0]), are zero.

In the present applications, both the sample size 𝒯 and the number of parameters ℒ are large: 𝒯≈(10^9) and ℒ≈(10^6). This makes direct estimation and forecasting difficult at best and infeasible at worst. The following material describes how use of polynomials and sparsity can help.

The first device we employ is due to Almon (1965), who suggests constraining the (Subscript[β, i]) to lie on polynomials in i. In the quadratic case, for example, (Subscript[β, i])=(Subscript[γ, 0])+(Subscript[γ, 1])i+(Subscript[γ, 2])(i^2). In matrix notation, this can be expressed as β=D γ, where γ=(({{(Subscript[γ, 0]), (Subscript[γ, 1]), (Subscript[γ, 2]))ᵀ}}) and D is the (ℒ×3)  design matrix:

D=(({{1, 1, 1}, {1, 2, 4}, {1, 3, 9}, { , ⋮,  }, {1, ℒ, (ℒ^2)}}))

Substituting into the regression equation gives

y=X β+e=X D γ+e

The first term on the right hand side can be grouped as X(D γ), which suggests a smoothing of the model parameters (the β coefficients). Grouping as (X D)γ suggests a transformation or smoothing of the data. This distinction, though, is purely a matter of perspective. The mathematics are the same.

A PDL can greatly reduce the dimensionality of the estimation. The transformed data matrix, X D is 𝒯×3, and the γ may be estimated directly as

(OverHat[γ])=(((DᵀXᵀX D))^(-1))DᵀXᵀy

This effectively reduces the size of the cross-product matrix: XᵀX is ℒ×ℒ, but DᵀXᵀX D is 3×3.

PDLs can also simplify the forecasting computations. Direct application of () requires ℒ multiplications and ℒ additions at each step. Using the PDL representation for the coefficients,

(Subscript[y, t])=((Subscript[γ, 0])+(Subscript[γ, 1])+(Subscript[γ, 2]))(Subscript[y, t-1])+((Subscript[γ, 0])+2(Subscript[γ, 1])+4(Subscript[γ, 2]))(Subscript[y, t-2])+⋯+((Subscript[γ, 0])+(Subscript[ℒγ, 1])+(ℒ^2)(Subscript[γ, 2]))(Subscript[y, t-ℒ])
=(Subscript[γ, 0])((Subscript[y, t-1])+(Subscript[y, t-2])+(Subscript[y, t-3])+⋯+(Subscript[y, t-ℒ]))+(Subscript[γ, 1])((Subscript[y, t-1])+2(Subscript[y, t-2])+3(Subscript[y, t-3])+⋯+(Subscript[ℒy, t-ℒ]))+(Subscript[γ, 2])((Subscript[y, t-1])+4(Subscript[y, t-2])+9(Subscript[y, t-3])+⋯+(ℒ^2)(Subscript[y, t-ℒ]))

Now define three state variables that will carry the quantities in the parentheses:

S(t,d)=(∑_(i=1)^ℒ)(i^d)(Subscript[y, t-i]), for d=0,1,2

To update S(t,d)->S(t+1,d), we first compute (Subscript[y, t])=(∑_(d=0)^2)(Subscript[γ, d])S(t,d). Then to take S(t,0)->S(t+1,0), note that (Subscript[y, t-ℒ]) leaves the sum and (Subscript[y, t]) enters:

S(t+1,0)=S(t,0)+(Subscript[y, t])-(Subscript[y, t-ℒ])

Similarly, by direct examination of the leaving and entering terms for higher orders:

S(t+1,1)=S(t,1)+S(t+1,0)-ℒ (Subscript[y, t-ℒ])
=S(t,1)+S(t,0)+(Subscript[y, t])-(ℒ+1)(Subscript[y, t-ℒ])
S(t+1,2)=S(t,2)+2S(t+1,1)-S(t+1,0)-(ℒ^2)(Subscript[y, t-ℒ])
=S(t,2)+2S(t,1)+S(t,0)+(Subscript[y, t])-(((1+ℒ))^2)(Subscript[y, t-ℒ])

It is still necessary to store all lagged values (Subscript[y, t]),..., (Subscript[y, t-ℒ]), but the update calculations are more efficient. (In the Matlab code, see the stateUpdate method in the polynom class.)

Additional notes on PDLs

The design matrix involves multiple polynomials covering different lagged intervals. These polynomials may be of different orders. Consider, for example, the design matrix

(({{1, 1, 0, 0, 0, 0}, {1, 2, 0, 0, 0, 0}, {1, 3, 0, 0, 0, 0}, {1, 4, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0}, {0, 0, 1, 0, 0, 0}, {0, 0, 1, 0, 0, 0}, {0, 0, 0, 1, 1, 1}, {0, 0, 0, 1, 2, 4}, {0, 0, 0, 1, 3, 9}}))

This specifies that: (Subscript[β, 1]),...,(Subscript[β, 4]) lie on a linear segment; (Subscript[β, 5])=(Subscript[β, 6])=(Subscript[β, 7]); and, (Subscript[β, 8]) through (Subscript[β, 10]) lie on a quadratic segment. Linear restrictions on the γ can be imposed to ensure that the segments join up. Splines can be used to give smoothness (continuity of derivatives, see Smith (1979)). Although PDLs are usually introduced using the direct polynomial forms used above, most estimation software packages use alternative representations that aim to control collinearity and errors of numerical precision (Cooper (1972)).

Sparsity

Application of PDL’s greatly reduces the size of the cross product matrix, from ℒ×ℒ to (in the quadratic case) 3×3. To arrive at that point in the direct way, though, we must compute X D where X is 𝒯×ℒ  and D is ℒ×3. This is still a formidable computation, so we take advantage of the sparsity of the data.  Specifically, if y and x are 𝒯-vectors of price changes, it may well be that the number of non-zero elements is on the order of (10^(-6))𝒯. This is sparsity in the usual sense, wherein a matrix is assumed to be primarily zeros, and only the non-zero elements (and their locations) are stored. Notationally, we write x={(Subscript[t, k]),(Subscript[v, k]) for k=1,...,K} where (Subscript[t, k]) and (Subscript[v, k]) denote the time and value of the (k^th) element. Implicitly, the K nonzero elements are embedded in a vector of size 𝒯, which has zeros (except at the (Subscript[t, k]) positions). Direct calculation of DᵀXᵀX D requires on the order of (𝒯^2)(ℒ^4) multiplications. By taking advantage of sparsity, this can be brought down to the order of (K^2).

Define ℐ(k) as a 𝒯×ℒ  matrix of zeros, with an ℒ×ℒ identity matrix located with its first element in row k and column 1. With this definition, the lagged data matrix may be written as

X=(∑_(k=1)^K)(Subscript[v, k])ℐ((Subscript[t, k])+1)

where the unit offset in (Subscript[t, k])+1 imposes a lag relative to the original series. It is understood that any terms with row index larger than 𝒯 are discarded.

For example, suppose that with 𝒯=12, there are three non-zero values of x: (Subscript[x, 1])=4, (Subscript[x, 3])=2, and (Subscript[x, 7])=5. Specified as a sparse matrix, t={1,3,7} and v={4,2,5} The lagged data matrix for ℒ=4 is X=4ℐ(2)+2ℐ(4)+5ℐ(8). Expressed in full form:

X=(({{0, 0, 0, 0}, {4, 0, 0, 0}, {0, 4, 0, 0}, {2, 0, 4, 0}, {0, 2, 0, 4}, {0, 0, 2, 0}, {0, 0, 0, 2}, {5, 0, 0, 0}, {0, 5, 0, 0}, {0, 0, 5, 0}, {0, 0, 0, 5}, {0, 0, 0, 0}}))=4×(({{0, 0, 0, 0}, {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}))+2×(({{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}))+5×(({{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}}))

It is initially easier to illustrate the calculations when the data are assumed to satisfy a non-overlap condition such that the (Subscript[t, k]) are sufficiently removed from one another that (Subscript[t, k])-(Subscript[t, k-1])≥ ℒ for k=2,...,K. In this case, the blocks of lagged variables are distinct. The example given above does not satisfy this condition. (Rows 4 and 5 contain overlapped lagged values.) A lagged data matrix that does satisfy the non-overlap condition can be constructed for nonzero values (Subscript[x, 1])=4, (Subscript[x, 4 ])=2, (Subscript[x, 7])=5, with 𝒯=10 and ℒ=3:

X=(({{0, 0, 0}, {4, 0, 0}, {0, 4, 0}, {0, 0, 4}, {2, 0, 0}, {0, 2, 0}, {0, 0, 2}, {5, 0, 0}, {0, 5, 0}, {0, 0, 5}}))=4×(({{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}}))+2×(({{0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}}))+5×(({{0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}}))

The cross product required for the OLS estimates is then:

DᵀXᵀX  D=Dᵀ((∑_k)(Subscript[v, k])ℐ((Subscript[t, k])+1)ᵀ)((∑_j)(Subscript[v, j])ℐ((Subscript[t, j])+1))D
=((∑_k)(Subscript[v, k])Dᵀℐ((Subscript[t, k])+1)ᵀ)((∑_j)(Subscript[v, j])ℐ((Subscript[t, j])+1) D)
=(∑_(k,j))(Subscript[v, j])(Subscript[v, k])Dᵀℐ((Subscript[t, k])+1)ᵀℐ((Subscript[t, j])+1) D
 =(∑_k)(Subscript[v, k]^2)Dᵀℐ((Subscript[t, k])+1)ᵀℐ((Subscript[t, k])+1) D

The last equality relies on the non-overlap condition. Since ℐ((Subscript[t, k])+1)ᵀℐ((Subscript[t, k])+1) is an identity matrix of order ℒ, the inner term of the last summation is equal to the product DᵀD, is a constant 3×3 matrix. Since we only need to compute it once, direct calculation is straightforward.  Of course, in this case, it might be pointed out that with the non-overlap condition, we could view the problem as consisting of K disjoint samples. Cumulating over these samples is the obvious thing to do, with no particular need to embed the data in the length-𝒯 time frame.

We now drop the non-overlap assumption. We also consider a more general bivariate cross product situation similar to those that arise in multivariate VARs. Suppose that both x and y are sparse time series. The lagged data matrix corresponding to x is X, a 𝒯×(Subscript[ℒ, x]) matrix, where ℒ denotes the number of lags on x. Similarly, the lagged data matrix for y is Y, a 𝒯×(Subscript[L, y]) matrix. The corresponding PDL design matrices are (Subscript[D, x]) and (Subscript[D, y]). (Subscript[D, x]) has (Subscript[ℒ, x]) rows, and (Subscript[D, y]) has (Subscript[ℒ, y]) rows. The orders of the polynomials (equal to the numbers of columns) may also differ.  The desired cross product matrix is (Subscript[D, x])ᵀXᵀY (Subscript[D, y]). The sparse representation of x involves time and value vectors (Subscript[t, k]) and (Subscript[v, k]) for k=1,...K≪𝒯; that of y has time and value vectors (Subscript[s, j]) and (Subscript[u, j]) for j=1,...,J≪𝒯. The lagged data matrices are represented as:

X=(∑_(k=1)^K)(Subscript[v, k])ℐ((Subscript[t, k])+1) and Y=(∑_(j=1)^J)(Subscript[u, j])𝒥((Subscript[s, j])+1)

where ℐ(t) is (as above) a 𝒯×(Subscript[ℒ, x]) matrix of zeros, with an identity matrix of order (Subscript[ℒ, x]) positioned with its leading element in the first column of row t, and 𝒥(s) is similarly defined as a 𝒯×(Subscript[ℒ, y]) matrix of zeros, with an identity matrix of order (Subscript[ℒ, y]) positioned at s. With these definitions,

(Subscript[D, x])ᵀXᵀY (Subscript[D, y])=(Subscript[D, x])ᵀ((∑_(k=1)^K)(Subscript[v, k])ℐ((Subscript[t, k])+1))ᵀ((∑_(j=1)^J)(Subscript[u, j])𝒥((Subscript[s, j])+1))(Subscript[D, y])

=(∑_(k=1)^K)(∑_(j=1)^J)(Subscript[v, k])(Subscript[u, j])(Subscript[D, x])ᵀℐ((Subscript[t, k])+1)ᵀ𝒥((Subscript[s, j])+1)(Subscript[D, y])

For convenience, define ℐ𝒥(t,s)≡ ℐ(t) ᵀ𝒥(s), a matrix of size (Subscript[ℒ, x])×(Subscript[ℒ, y]), and with elements are to zero or one. ℐ𝒥(t,s) indicates the overlap between the two lagged structures. Suppose, for example, that with 𝒯=10, (Subscript[ℒ, x])=3, (Subscript[ℒ, y])=4, (Subscript[t, 1])=2, and (Subscript[s, 1])=3. Then

ℐ(2)=(({{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}, {0, 0, 0}})),𝒥(3)=(({{0, 0, 0, 0}, {0, 0, 0, 0}, {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}})), and ℐ𝒥(2,3)=(({{0, 0, 0, 0}, {1, 0, 0, 0}, {0, 1, 0, 0}}))

These representations can greatly simplify the cross product calculation. Initially suppose that each of the PDL design matrices contains only a constant term: (Subscript[D, x])=(({{1, 1, 1}}))ᵀ and (Subscript[D, y])=(({{1, 1, 1, 1}}))ᵀ. Then the contribution to the cross product is

(Subscript[v, 1])(Subscript[u, 1])(({{1, 1, 1}}))(({{0, 0, 0, 0}, {1, 0, 0, 0}, {0, 1, 0, 0}}))(({{1}, {1}, {1}, {1}}))=(Subscript[v, 1])(Subscript[u, 1])×2

If the PDL design matrices contain constant and linear terms,

(Subscript[D, x])ᵀ=(({{1, 1, 1}, {1, 2, 3}}))and (Subscript[D, y])ᵀ=(({{1, 1, 1, 1}, {1, 2, 3, 4}})),

and the contribution to the cross product is:

(Subscript[v, 1])(Subscript[u, 1])(({{1, 1, 1}, {1, 2, 3}}))(({{0, 0, 0, 0}, {1, 0, 0, 0}, {0, 1, 0, 0}}))(({{1, 1}, {1, 2}, {1, 3}, {1, 4}}))=(Subscript[v, 1])(Subscript[u, 1])(({{2, 3}, {5, 8}}))

With quadratic PDL design matrices the cross product contribution is

(Subscript[v, 1])(Subscript[u, 1])(({{1, 1, 1}, {1, 2, 3}, {1, 4, 9}}))(({{0, 0, 0, 0}, {1, 0, 0, 0}, {0, 1, 0, 0}}))(({{1, 1, 1}, {1, 2, 4}, {1, 3, 9}, {1, 4, 16}}))=(Subscript[v, 1])(Subscript[u, 1])(({{2, 3, 5}, {5, 8, 14}, {13, 22, 40}}))

In general, the entries in the last matrix are double partial sums of offset finite geometric series, for which compact closed-form expressions are readily available. (It is not necessary to construct and multiply the matrices on the left hand sides of these expressions. See the Mathematica section on lagged polynomial cross products of sparse vectors.) Taking advantage of these closed-form expressions, define f(t,s,(Subscript[ℒ, x]),(Subscript[ℒ, y]),(Subscript[d, x]),(Subscript[d, y]))=(Subscript[D, x])ᵀ ℐ(t)𝒥(s)(Subscript[D, y]), where (Subscript[D, x]) is the PDL design matrix for x, with maximum degree (Subscript[d, x]) and lag length (Subscript[ℒ, x]). (Subscript[D, y]) is defined similarly.

(Subscript[D, x])ᵀXᵀY (Subscript[D, y])=(∑_(k=1)^K)(∑_(j=1)^J)(Subscript[v, k])(Subscript[u, j])f((Subscript[t, k])+1,(Subscript[s, j])+1,(Subscript[ℒ, x]),(Subscript[ℒ, y]),(Subscript[d, x]),(Subscript[d, y]))

Direct computation of the left-hand side involves something on the order of (Subscript[(ℒ^2), x])×(𝒯^2)×(Subscript[(ℒ^2), y]) operations. Evaluation of finvolves (at most) several hundred operations. The double sum is only over the nonzero values of x and y. If we sort the entries in the sparse matrices by time (of the nonzero elements), we can realize greater efficiencies. For a given k, the values of j that make nonzero contributions to the crossproduct (that is, ℐ𝒥((Subscript[t, k])+1,(Subscript[s, j])+1)≠0) can be determined before the start of the inner loop.

References

Almon, Shirley, 1965, The Distributed Lag Between Capital Appropriations and Expenditures, Econometrica 33, 178-196.

Cooper, J. Phillip, 1972, Two Approaches to Polynomial Distributed Lags Estimation: An Expository Note and Comment, The American Statistician 26, 32-35.

Smith, Patricia L., 1979, Splines as a useful and convenient statistical tool, American Statistician 33, 57-62.

Mathematica Code

Polynomial update verifications

This section may be evaluated.

Clear[y,S,t,n,d];

The PDL state variables are as follows. t is the current time index; n is the max lag; d is the degree.

S[t_:t,n_:n,d_:0]:=Sum[(Subscript[y, t-i])(i^d),{i,1,n}]

{S[t,n,#]&/@{0,1,2}}//TableForm

Verify the degree 0 update expression for various values of n:

Table[Evaluate[S[t+1,n,0]==(S[t,n,0]+(Subscript[y, t])-(Subscript[y, t-n]))],{n,5,20,5}]//Simplify

Verify the degree 1 update expression for various values of n.

Table[S[t+1,n,1]==(S[t,n,1]+S[t+1,n,0]-n (Subscript[y, t-n])),{n,5,20,5}]

Alternatively, the expression used in the Matlab code is:

(S[t,n,1]+(S[t,n,0]+(Subscript[y, t])-(Subscript[y, t-n]))-n (Subscript[y, t-n]))//Simplify

Table[S[t+1,n,1]==(S[t,n,1]+(S[t,n,0]+(Subscript[y, t]))-(n+1) (Subscript[y, t-n])),{n,5,20,5}]

Verify the degree 2 update expression for various values of n.

Table[S[t+1,n,2]==S[t,n,2]+2S[t+1,n,1]-S[t+1,n,0]-(n^2)(Subscript[y, t-n]),{n,5,20,5}]//Simplify

Alternatively, the expression used in the Matlab code is:

S[t,n,2]+2(S[t,n,1]+(S[t,n,0]+(Subscript[y, t]))-(n+1) (Subscript[y, t-n]))-(S[t,n,0]+(Subscript[y, t])-(Subscript[y, t-n]))-(n^2)(Subscript[y, t-n])//Simplify

Table[S[t+1,n,2]==S[t,n,2]+2(S[t,n,1]+(S[t,n,0]+(Subscript[y, t]))-(n+1) (Subscript[y, t-n]))-(S[t,n,0]+(Subscript[y, t])-(Subscript[y, t-n]))-(n^2)(Subscript[y, t-n]),{n,5,20,5}]//Simplify

Lagged polynomial cross products of sparse vectors

Matrix construction

Clear[X,Y,𝒟,tSum];

𝒟 denotes the coefficient (design) matrix for polynomials. (“D” is reserved in Mathematica for derivatives.)

𝒟[n_,d_]:=Table[(i^dd),{dd,0,d},{i,1,n}]\[Transpose];
𝒟[4,2]//MatrixForm

X[i,t] is a matrix with 𝒯 rows, n columns and an identity matrix starting in row i. It represents a lag matrix, where a variable that appears in position i,1 also appears in position i-1,2, position i-2,3, ... position i-n+1,n. (The Mathematica Block construct is used to localize the assignment 𝒯=8.)

X[i_,n_]:=Module[{x},
x=Table[0,{𝒯},{n}];
If[i+n-1>𝒯,Print["X[i,n]: i+n-1>𝒯 (i=",i,"; n=",n,"; 𝒯=",𝒯,")"];Return[Null]];
x[[i;;i+n-1]]=IdentityMatrix[n];
x];
Block[{𝒯=8},X[3,4]]//MatrixForm

Define Y similarly:

Y[i_,n_]:=X[i,n]

and

Ones:=Table[{1},{𝒯}];
Block[{𝒯=8},Ones]//MatrixForm

Matrices used in text examples:

Block[{𝒯=12,ℒ=4},x1=X[2,ℒ];x2=X[4,ℒ];x3=X[8,ℒ];r=4*x1+2*x2+5*x3;MatrixForm/@{x1,x2,x3,r}]

Block[{𝒯=10,ℒ=3},x1=X[2,ℒ];x2=X[5,ℒ];x3=X[8,ℒ];r=4*x1+2*x2+5*x3;MatrixForm/@{r,x1,x2,x3}]

Block[{𝒯=10},x=X[2,3];y=X[3,4];xty=x\[Transpose].y;MatrixForm/@{x,y,xty}]

MatrixForm/@{𝒟[3,1],𝒟[4,1],𝒟[3,1]\[Transpose].(({{0, 0, 0, 0}, {1, 0, 0, 0}, {0, 1, 0, 0}})).𝒟[4,1]}

MatrixForm/@{𝒟[3,2],𝒟[4,2],𝒟[3,2]\[Transpose].(({{0, 0, 0, 0}, {1, 0, 0, 0}, {0, 1, 0, 0}})).𝒟[4,2]}

Continuation

The result matrix is (Subscript[𝒟, x])ᵀ.Xᵀ.Y.(Subscript[𝒟, y]). The following depicts the pieces and final result in the case of fixed values for 𝒯, ℒx, and ℒy; for X[5,ℒx] and Y[j,ℒy] where j goes from 2 to 9.

Block[{𝒯=12,ℒx=4,ℒy=3},
Table[{j,
X[5,ℒx]//MatrixForm,
Y[j,ℒy]//MatrixForm,
X[5,ℒx]\[Transpose].Y[j,ℒy]//MatrixForm,
𝒟[ℒx,2]\[Transpose].X[5,ℒx]\[Transpose].X[j,ℒy].𝒟[ℒy,2]//MatrixForm},{j,2,9}]//Transpose//TableForm]//Style[#,FontSize->9]&

Depending on the relative offset between the X and Y matrices, the XᵀY product is a diagonal band of ones that moves from top/right to bottom/left.

Here is an alternative calculation. i indicates the location of the identity matrix in X, and j indicates the location for Y.

CP0[i_,j_,ℒx_,ℒy_]:=Module[{k,ka,kb,z},
z=Table[0,{3},{3}];
k=j-i; (*offset*)
ka=Max[1,k+1];
kb=Min[ℒx,ℒy+k];
If[ka>kb,Return[z]];
Table[Sum[(ii^d1)(((ii-k))^d2),{ii,ka,kb}],{d1,0,2},{d2,0,2}]];
Table[{j,CP0[5,j,4,3]//MatrixForm},{j,2,9}]//Transpose//TableForm//Style[#,FontSize->10]&

In CP0 the geometric sums in the matrix given by

Table[Sum[(ii^d1)(((ii-k))^d2),{ii,ka,kb}],{d1,0,2},{d2,0,2}]

are evaluated in full form (with the number of iterations approximately equal to Min[ℒx, ℒy]). It is very useful to have closed-form representations for these sums. Below is a grid with rows indexed by d1 and columns indexed by d2.

Clear[i,k,ka,kb,ia,ib];

t=Table[Sum[(i^d1)(((i-k))^d2),{i,ia,ib}],{d1,0,2},{d2,0,2}];

Grid[t,Frame->All,ItemStyle->{FontSize->12},ItemSize->13,Alignment->{Center,Center}]

(In the Matlab code, these calculations are implemented in the pSum method of the polynom class.)

Modifications and extensions used in Matlab code

Adjusting the lagged data matrix to conform to [firstValid, lastValid]

In the Matlab programing, the data vectors are implicitly defined on i=1,...,𝒯.  The lagged data matrix X is assumed to have 𝒯 rows. (In the VAR applications, 𝒯 is the number of time units per day, 3600×24=86,400 if the time units are seconds, and  86,400,000 if the time units are milliseconds.) Due to lags and incomplete data, though, not all rows of X are valid. The valid range is indicated as 0<firstValid<lastValid≤ 𝒯. All routines are adjusted to confine calculations to these rows. This section uses the following conventions.

i is the original position in the sparse data vector corresponding to x (i.e., at lag zero).

a is an offset relative to i, where the PDL starts.

Similarly, j and b are the original position for y and the offset for the y PDL.

vTrim removes all rows from X that are outside of [firstValid, lastValid].

vTrim[X_]:=Module[{d},
d=Dimensions[X];
If[Not[MatrixQ[X]] || lastValid>d[[1]] || firstValid<1||firstValid>lastValid,
Print["vTrim error."];Return[Null]];
X[[firstValid;;lastValid]]];
Block[{𝒯=10,firstValid=2,lastValid=8},MatrixForm/@{X[3,4],vTrim[X[3,4]]}]

Here is a cross product calculation that uses vTrim to adjust the range.

Block[{𝒯=20,firstValid=8,lastValid=10,j=10,a=0,b=0,ℒx=4,ℒy=3},
Table[{
i,
(*X[i+a,ℒx]//MatrixForm,
Y[j+b,ℒy]//MatrixForm,
vTrim[X[i+a,ℒx]]//MatrixForm,
vTrim[Y[j+b,ℒy]]//MatrixForm,*)
X[i+a,ℒx]\[Transpose].Y[j+b,ℒy]//MatrixForm,
vTrim[ X[i+a,ℒx] ]\[Transpose].vTrim[ Y[j+b,ℒy] ]//MatrixForm,
𝒟[ℒx,2]\[Transpose].vTrim[X[i+a,ℒx]]\[Transpose].vTrim[Y[j+b,ℒy]].𝒟[ℒy,2]//MatrixForm},{i,6,13}]//Transpose//TableForm]//Style[#,FontSize->9]&

The code used in Matlab adjusts the starting and ending values of the loop to conform to firstValid and lastValid:

CP1[i_,j_,a_,b_,n_,m_,dx_,dy_]:=Module[{k,ka,kb,z},
z=Table[0,{dx+1},{dy+1}];
k=(j+b)-(i+a); (*offset*)
ka=Max[1,k+1,firstValid-(i+a)+1];
kb=Min[n,m+k,lastValid-(j+b)+1+k];
If[ka>kb,Return[z]];
Table[Sum[(ii^d1)(((ii-k))^d2),{ii,ka,kb}],{d1,0,dx},{d2,0,dy}]];
Block[{𝒯=20,firstValid=8,lastValid=10,j=10,a=0,b=0,ℒx=4,ℒy=3},
Table[{i,CP1[i,j,a,b,ℒx,ℒy,2,2]//MatrixForm},{i,j-4,j+3}]//Transpose//TableForm]//Style[#,FontSize->9]&

Verify equivalence of calculation methods for a range of values ...

kStep=0;
Block[{𝒯=30},
Do[
Block[{firstValid=kfv,lastValid=klv,i=ki,j=kj,a=kka,b=kkb,dx=2,dy=1,n=4,m=1},
kStep = kStep+1;
If[Mod[kStep,2500]==0 ,Print[kStep," ",kfv, " ", klv," ",ki," ",kj," ",kka]];
If[kStep>15000,Break[]];
z1=CP1[i,j,a,b,n,m,dx,dy];
z2=𝒟[n,dx]\[Transpose].vTrim[X[i+a,n]]\[Transpose].vTrim[X[j+b,m]].𝒟[m,dy];
(*z2=123;*)
If[z1===z2,Continue[],Print["Whoops"]Print[z1//MatrixForm,z2//MatrixForm];Break[];];],
{kfv,5,15,5},{klv,kfv+10,𝒯-5,5},{ki,Max[1,kfv-5],𝒯-8},{kj,Max[1,kfv-5],𝒯-8,2},{kka,0,2},{kkb,0,2}]];

Non-sparse data series

The sparse data vectors discussed to this point have all been sparse (like price changes). A price discovery VECM will also often include a constant term and/or a cointegration vector. These require special handling.

Constant terms

Many models will include constant terms. In the data matrix, this is usually represented as a column of ones. The constant crossed with itself is simply the number of observations, N. We might have N=𝒯, the length of the implicit data matrix. More generally, though, N=lastValid-firstValid+1, where firstValid and lastValid are (as described earlier) pointers in X.

The cross product of the constant vector and a sparse data vector is the sum of the nonzero values.

Letting (Subscript[ι, l]) denote a column vector of l ones, the cross product of the constant vector and a quadratic PDL on a sparse data vector is:

(Subscript[D, x])ᵀXᵀ(Subscript[ι, 𝒯])=(Subscript[D, x])ᵀ((∑_(k=1)^K)(Subscript[v, k])ℐ((Subscript[t, k])+1))ᵀ(Subscript[ι, 𝒯])=(∑_(k=1)^K)(Subscript[v, k])(Subscript[D, x])ᵀℐ((Subscript[t, k])+1)ᵀ(Subscript[ι, 𝒯])=(∑_(k=1)^K)(Subscript[v, k])(Subscript[D, x])ᵀ(Subscript[ι, ℒ ])

where ℒ is the length of the PDL. The (Subscript[D, x])ᵀ(Subscript[ι, ℒ ]) term is the sum (across columns) of the elements of (Subscript[𝒟, x])ᵀ.

Block[{𝒯=10,ℒx=4},{MatrixForm[#]&/@{
X[5,ℒx],
X[5,ℒx]\[Transpose].Ones,
𝒟[ℒx,2]\[Transpose],
𝒟[ℒx,2]\[Transpose].X[5,ℒx]\[Transpose].Ones}}]//TableForm

Sparse levels

Cointegrating vectors are not sparse in the usual sense (of having a small number of nonzero values). The difference between two price series (Subscript[Δ, t])=(Subscript[p, 1t])-(Subscript[p, 2t]) is not sparse, even if the price changes of component prices are sparse. (Subscript[Δ, t]) is sparse in a different sense, sometimes called piecewise sparse or level sparse. Its entries consist of constant levels that persist in time until one of component prices changes. The cross product calculations are adjusted accordingly.

Matlab notes

The Matlab code has two sparse matrix classes. spd matrices, used for price differences, are sparse in the usual sense. spl matrices are level sparse.

A VECM analysis involves four types of data series:

The Crossproduct class has static methods (functions) to compute crossproducts involving all pairwise combinations of these types.