# Why is EstimatedProcess so slow for a vector ARProcess?

I'm trying to fit some 4 x N data to a vector auto-regressive model. However, EstimatedProcess is very slow for this. For a 4-th order vector AR process, it takes more than 10 seconds (for N approx. equal to 2500). When I try a 5th order AR process, it takes 185 seconds. I noticed that, by default, Mathematica apparently uses the "MethodOfMoments" method for EstimatedProcess applied to ARProcess. So, I coded up my own function for fitting a vector AR process, using the Method of Moments:

VARProcessMOM[data_List, order_Integer] :=
Module[{dataPrec, i, j, k, n, mn, lag, dataD, cov, covVec, covMat, covMatPar, param, covar, useBiasedEst = True, precision = 50},
$MinPrecision = precision; (* 'data' dimensions must be k x n, where: k = number of arrays, n = length of arrays *) {k, n} = Dimensions[data]; dataPrec = SetPrecision[data, precision]; mn = Mean[#] & /@ dataPrec; dataD = Transpose[dataPrec - mn]; If[useBiasedEst, (* Mathematica EstimatedProcess[data, ARProcess[...], Method->”MethodOfMoments] apparently uses n in the denominator of the autocovariance calculation. (See 'CovarianceFunction' definition in Mathematica.) This produces a biased estimate. Nevertheless, this produces results nearly identical to Mathematica. *) cov = Table[ Sum[KroneckerProduct[dataD[[i]], dataD[[i + lag]]], {i, 1, n - lag}]/n, {lag, 0, order} ], (* Using (n-lag) in the denominator gives an unbiased estimate of the autocovariance *) cov = Table[ Sum[KroneckerProduct[dataD[[i]], dataD[[i + lag]]], {i, 1, n - lag}]/(n - lag), {lag, 0, order} ] ]; covVec = ArrayFlatten[Transpose[{cov[[2 ;;]]}]]; covMat = Table[ If[j > i, Transpose[cov[[j - i + 1]]], cov[[i - j + 1]]], {i, 1, order}, {j, 1, order} ]; covMat = ArrayFlatten[covMat]; covMatPar = Inverse[Transpose[covMat].covMat].Transpose[covMat].covVec; param = Transpose[#] & /@ Partition[covMatPar, k]; covar = cov[[1]] - Sum[param[[i]].cov[[i + 1]], {i, 1, order}];$MinPrecision = MachinePrecision;
param = SetPrecision[param, MachinePrecision];
covar = SetPrecision[covar, MachinePrecision];
{param, covar}
]


This function produces results that are essentially identical to EstimatedProcess/ARProcess (when "useBiasedEst" is set to True). For scalar input data, you can just enclose the data in brackets, {data}, so that it becomes a 1 x N matrix. The input "order" is just the number of desired autoregressive terms (lags) to be included.

For example, the output produced by EstimatedProcess for a 2-vector AR process with two lags:

r1 = {{.9, .3}, {-.2, .8}};
r2 = {{.1, -.02}, {.05, .05}};
r1 // MatrixForm
r2 // MatrixForm
cov = {{.3, .02}, {.02, .5}};
cov // MatrixForm
RandomSeed[9];
data = RandomFunction[ARProcess[{r1, r2}, cov], {1, 10^2}];
es1 = EstimatedProcess[data, ARProcess[Array[aa, {2, 2, 2}], Array[ee, {2, 2}]]];
es1[[1]]
es1[[2]]


is:

{{{0.918574, 0.321707}, {-0.253664, 0.641996}}, {{0.0466231, 0.0096257}, {0.0951791, 0.222094}}}

{{0.413422, 0.0454421}, {0.0454421, 0.515956}}


This is identical to the output produced by my VARProcessMOM function (when useBiasedEst is set to True):

es2 = VARProcessMOM[Transpose[Normal[data][[1, All, 2]]], 2];
es2[[1]]
es2[[2]]


However, when I take the differences of these outputs, there is a small difference in the numbers between EstimatedProcess[data, ARProcess[...]]] and VARPRocessMOM, on the order of 10^-11 or smaller. You will notice that I increased the precision of the calculations within the program. Previous to doing this, the differences were slightly larger. Also, for large data sets and large lags (4 or greater), I was getting the "badly conditioned matrix" error when taking the matrix inverse. Increasing the precision of calculations solved this problem. The small remaining differences between the outputs are not producing any discernible difference in my application: forecasting 1-step-ahead prices for a weakly stationary portfolio of exchange-traded funds. I'm getting good results. I just wanted to speed up the calculations.

The big difference between the two functions is in the run time. For the 4 x 2500 data set, it took 10.7 seconds for the EstimateProcess function (using order = 4). For my VARProcessMOM function it only took 0.18 seconds to get the same result. I can now experiment with using more lags, because it doesn't slow down the calculations significantly.

Why is the Mathematica EstimatedProcess function so slow? Are they using more than 50 digits of precision? That seems unnecessary. I don't see any change in the output of my function after about 40 digits. I suppose for really large data sets, or large numbers of lags, perhaps it might make a difference. However, the EstimatedProcess function is so slow, you can't process large data sets with many lags anyway.