function [W, C, D, R, S, information, rho, elapsed, iterations, primalResidual, dualResidual, logLikelihood, stepsizes, iterationTimes, logLikelihoods, fixedIterations, handles, gapHandle, dLikelihood] = ...
  bst_vgarch(model, epsilon, nTrials, R, S, flagELM, method, qualifier, maxIterations, initial, rho, VERBOSE, ax1, lineColor, hLikelihood, handlePrefix)
% BST_VGARCH        Fit a vector GARCH model to the signal with chosen order
%                   in conditional heteroscedasticity and
%                   autoregression-in-variance.
% Inputs:
%   model           - model used for vector GARCH modelling
%                     'BEKK', 'bekk': use symmetric, constrained BEKK model
%                      'vec',  'vec': use more general, simpler vec model
%                     ^-- by default, the 'vec' model is used
%                     [Z: string]
%   e               - vector signal that is autoregressive in variance
%                     [E: M x NT double]
%   nTrials         - number of trials in the signal
%                     [T: scalar positive integer]
%   R               - order of autoregression-in-variance
%                     [R: scalar nonnegative integer]
%   S               - Order of heteroscedasticity-in-variance
%                     [S: scalar nonnegative integer]
%                     ^-- currently, only S = 0 works
%   flagELM         - if true, use Engle's Lagrange multiplier statistic to
%                     optimize R in the ARCH model
%   method          - method used to estimate parameters
%                     1: quadratic programming
%                     2: preconditioned conjugate gradient descent
%                     3: alternating-direction method of multipliers
%                     other: error
%                     [C: positive integer]
%   qualifier       - desired modified for the method chosen above in C
%                     For quadratic programming,
%                     'off': use numerical gradient
%                      'on': use analytic gradient
%                     For preconditioned conjugate gradient descent,
%                     'off': use numerical gradient
%                      'on': use analytic gradient
%                     For alternating direction method of multipliers,
%                     the qualifier is not used.
%                     [Q: string]
%   maxIterations   - maximum number of iterations for any of the methods
%                     [I: scalar positive integer]
%   initial         - initial vector of squares of values, stacked into a column
%                     For BEKK, (remember here C/D are the small matrices)
%                       [lower triangle of W; C; D]
%                     For vec, (remember here C/D are the big matrices)
%                       [W; C; D]
%                     If a scalar, then
%                        1: use MFEToolbox (I don't know what it does)
%                     BEKK: use scalar regression for each signal
%                      vec: use multivariate regression on all pairs of signals
%                     [X0: vector nonnegative double]
%   rho             - augmented Lagrangian parameter for ADMM
%                     This parameter is not used for the other methods (C = 1,2)
%                     [P: scalar positive double]
%   VERBOSE         - flag to show iterative results in optimization methods
%                     [V: boolean]
%   ax1             - likelihood function plot handle for showing traversal
%                     If nonempty, we plot the traversal for each iteration onto
%                     the axis handle
%                     [AX: scalar positive double corresponding to axis handle]
%   lineColor       - color of traversal line
%                     [LC: 3 x 1 double w/ each element in [0,1]]
%   hLikelihood     - handle to figure of plotted log-likelihood, used to save
%                     each traversal in a frame. Keep empty to save nothing.
%                     [HL: scalar positive double corresponding to axis handle]
%   handlePrefix    - prefix to attach to file when saving image, if the
%                     likelihood handle is given
%                     [H: string]
% Outputs:
%   W               - homoscedastic covariance term in GARCH model
%                     [W: M x M or M^2 x 1 positive double]
%   C               - conditional heteroscedasticity coefficients in GARCH model
%                     [C: M x M x R or M^2 x M^2 x R positive double]
%   D               - autoregression-in-variance coefficients in GARCH model
%                     [D: M x M x S or M^2 x M^2 x S positive double]
%   R               - lag used for conditional heteroscedasticity
%   S               - lag used for variance autoregression
%   information     - information matrix from log-likelihood for parameter
%                     estimates
%   rho             - aDMM parameter that provided convergence, if ADMM was chosen
%   elapsed         - time to compute coefficients, in seconds
%   iterations      - bumber of iterations
%   primalResidual  - primal residual as a fcn of iterations
%   dualResidual    - dual residual as a fcn of iterations
%   logLikelihood   - average log-likelihood for the estimated parameters
%   stepsizes       - change in parameters at each iteration of method chosen
%   iterationTimes  - elapsed time for each iteration of method chosen
%   logLikelihoods  - average log-likelihood for each iteration of method chosen
%   fixedIterations - number of fixed-point runs in each ADMM iteration
%   handles         - handles to all line segments in the traversal, if plotted
%   gapHandle       - handle to plot of duality gap
%   dLikelihood     - change in likelihood at each iteration

% @=============================================================================
% This software is part of the Brainstorm software:
% http://neuroimage.usc.edu/brainstorm
% 
% Copyright (c)2000-2013 Brainstorm by the University of Southern California
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html.
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Syed Ashrafulla, 2012

% This code uses mmx, a software package that allows many matrix multiplications
% at the same time. To do this, you need to download build_mmx from
% http://www.mathworks.com/matlabcentral/fileexchange/37515-mmx-multithreaded-matrix-operations-on-n-d-matrices
% and then go to that folder and edit build_mmx.m to modify one line and add two lines:
% -- on line 80, change the iteration to go from 2 to 3
% ---- do not do this if you have the Intel MKL installed. You won't have this
% ---- installed unless you know it immediately because you are a developer.
% -- on line 196, add "syed = pwd"
% -- on line 214, add "cd(syed)"
% Next, you need to update /Applications/MATLABR2013a.app/bin/mexopts.sh. Use
% vim and the command
% :%s/10.7/10.8
% and you will get 5 substitutions. Write it out using :wq! to force overwriting
% and then go back into MatLab and hit mex -setup. Pick 1 to get your edits onto
% your use of MatLab.
% Finally, build_mmx will work with 1 stupid type conversion warning!

% Some notation help here, it follows the math on the website:
% M = dimension of signal
% N = number of valid timepoints in the sense that a log-likelihood can be calculated (so at each trial, starting at max(R,S)+1 and going to last sample)
% e = vectorized local covariation, i.e. epsilon_{n-1} epsilon_{n-1}^T, which is the "data" used in the modelling (in bst_garch, this is e.^2)
% eMatrix = unvectorized version of e, which is used in the ADMM approach; since it is outside of iteration, it is of no cost to the other methods

totalStart = tic;

% dimension of innovation
M = size(epsilon,1); N = size(epsilon,2); tolerance = 1e-6;
nOutputs = nargout;

% "squared" innovations
eMatrix = zeros(M, M, N); % this is the unvectorized "e", which I use in the ADMM part (the R_n)
for n = 1:N
  eMatrix(:,:,n) = epsilon(:, n) * epsilon(:, n)';
end
e = reshape(eMatrix, M*M, N); % vectorize to get the "e" in the notes

% interpolation for time indices below 0
weights = .06 * .94.^(0:(ceil(sqrt(N))-1)); weights = weights/sum(weights);
backCast = sum(bst_bsxfun(@times, weights, e(:,1:ceil(sqrt(N)))), 2);
e = reshape([repmat(backCast, [1 max(R,S) nTrials]) reshape(e, [M*M N/nTrials nTrials])], M*M, max(R,S)*nTrials+N);
eMatrix = reshape(e, M, M, max(R,S)*nTrials+N);

% convert full-vectorized squared innovations to half-vectorized square innovations
e = e(tril(true(M)), :);

% timepoints
flagMultiple = (nTrials > 1);
if flagMultiple
  n = bst_trial_idx(max(R,S)+1 : size(e,2)/nTrials, size(e,2)/nTrials, nTrials)';
else % Send empty array so bst_mvar_transfer knows not to call a fcn
  n = max(R,S)+1 : size(e,2);
end

% assignment from lower-triangular half-vectorization to full-vectorization and back
toFull = zeros(M); toFull(tril(true(M))) = 1:(M*(M+1)/2); toFull = tril(toFull, -1)' + toFull; toFull = toFull(:);
toHalf = repmat(tril(true(M)), [1 1 N]);

% model choice from string to flag: BEKK is true
model = strcmpi(model, 'BEKK');

%% Optimal order for conditional heteroscedasticity modeling
% we adjust R based on Engle's Lagrange multiplier statistic
% = N * (1 - SS_(R lags) / SS_(R-1 lags)) which is chi-square with two degrees of freedom for our bivariate systems

if flagELM
  
  % pre-allocate the statistic
  engle = zeros(1, R+1);
  
  % lag-o testing (so whether removing the mean is useful)
  base = squeeze(epsilon(:, n)).^2; % residual without any modelling
  SS0 = sum(sum(( base ).^2, 2)); % H_0: no regression whatsoever
  coefficients = ones(1, size(base,2)) * mean(base(:)); % coefficients only remove the mean
  SS1 = sum(sum(( base - (base / coefficients) * coefficients ).^2)); % H_1: regress out a constant (the "mean" single-sample second moment) from signals
  engle(1) = size(base,2) * (1 - SS1/SS0); % Engle's Lagrange multiplier statistic: chi-square with 2 degrees of freedom

  % regression testing on the sum of squares for the remaining lags iteratively (based on the previous lag)
  for r = 1:R
    SS0 = SS1; % H_0: regression with R-1 lags
    coefficients(end+(1:2),:) = squeeze(epsilon(:, n-r)).^2; % add the Rth lag for both signals
    SS1 = sum(sum(( base - (base / coefficients) * coefficients ).^2)); % H_1: regression with R lags
    engle(r+1) = size(base,2) * (1 - SS1/SS0); % Engle's Lagrange multiplier statistic: chi-square with 2 degrees of freedom
  end
  
  % convert to maximum lag in ARCH modeling
  R = max(0, sum(engle > 5.99146454710798082032852, 3)-1); % all lags with statistic > chi2inv(1-0.05,2) are admitted, with min lag zero
  
end

%% For no conditional heteroscedasticity model, pre-send the results
if R == 0
  
  % conditional homoscedasticity so C = D = 0
  W = bst_correlation(epsilon, [], struct('normalize', false, 'nTrials', 1, 'maxDelay', 0, 'nDelay', 1, 'flagStatistics', false));
  if model % for BEKK, we have
    C = zeros(M, M, R); % the conditional heteroscedasticity coefficients (C * e * e^T * C)
    D = zeros(M, M, S); % the autoregression-in-variance coefficients (D * h * h^T * D^T)
    information = eye( M*(M+1)/2 + R*M*M + S*M*M);
  else % for vec, we have
    C = zeros(M*(M+1)/2, R*M*(M+1)/2); % the conditional heteroscedasticity coefficients (C * vech(e))
    D = zeros(M*(M+1)/2, S*M*(M+1)/2); % the autoregression-in-variance coefficients (D * vech(h))
    information = eye( M*(M+1)/2 + R*M*(M+1)/2*M*(M+1)/2 + S*M*(M+1)/2*M*(M+1)/2);
  end
  
  % runtime
  elapsed = toc(totalStart);
  if VERBOSE
    fprintf(1, 'Elapsed time is %0.6f seconds.\n', elapsed);
  end

  % all the other outputs
  iterations = 0;
  primalResidual = 0;
  dualResidual = 0;
  logLikelihood = Inf;
  stepsizes = [];
  iterationTimes = [];
  logLikelihoods = [];
  fixedIterations = [];
  handles = [];
  gapHandle = [];
  dLikelihood = [];
  
  % hard return
  return;
  
end

%% Initialization
if size(initial,1) > 1 % set the parameters to those in the initialization

  % assign for ADMM method
  W = initial(1:M*(M+1)/2); % the homoscedastic covariance's lower half-triangle
  if model % for BEKK, we have
    C = reshape(initial(M*(M+1)/2 + (1:R*M*M)), M, M, R); % the conditional heteroscedasticity coefficients (C * e * e^T * C)
    D = reshape(initial(M*(M+1)/2 + R*M*M + (1:S*M*M)), M, M, S); % the autoregression-in-variance coefficients (D * h * h^T * D^T)
  else % for vec, we have
    C = reshape(initial(M*(M+1)/2 + (1:R*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, R*M*(M+1)/2); % the conditional heteroscedasticity coefficients (C * vech(e))
    D = reshape(initial(M*(M+1)/2 + R*M*(M+1)/2*M*(M+1)/2 + (1:S*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, S*M*(M+1)/2); % the autoregression-in-variance coefficients (D * vech(h))
  end

else % try to find an initialization that's decent
  
  % first try least-squares
%   [W, C, D] = bst_vgarch_initialize(model, e, n, M, R, S, N, toFull, toHalf, 3, '', 999, tolerance, false);
%   if isinf(bst_vgarch_likelihood(model, [W; C(:); D(:)], M, R, S, e, n, toFull, toHalf, false)) % if least-squares fails to get a feasible point, ...
%     [W, C, D] = bst_vgarch_initialize(model, e, n, M, R, S, N, toFull, toHalf, 2, 'off', 999, tolerance, false); % ... then do scalar variance targeting
%   end
  
  % forget least squares and use scalar variance targeting
  [W, C, D] = bst_vgarch_initialize(model, e, n, M, R, S, N, toFull, toHalf, 2, 'off', 999, tolerance, false);
  
end

% save the first sum-of-squares for plotting purposes
if ~isempty(ax1)
  oldW2 = sum(W(:).^2); oldC2 = sum(C(:).^2); % oldD2 = sum(D(:).^2);
end
% this is done again in the iteration loop, but now I've forced the scope of the variable to the parent function.
% that scope increase will mean that in the next iteration I now have the old sum-of-squares saved here.

%% Testing
% nParameters = M*(M+1)/2 + R*M*(M+1)/2*M*(M+1)/2 + S*M*(M+1)/2*M*(M+1)/2;
% x = [W(:); C(:); D(:)] + [rand(M*(M+1)/2,1) * 0.2 - 0.1; rand(R*M*(M+1)/2*M*(M+1)/2, 1)*0.04 - 0.02];
% ground = zeros(nParameters, 1); [~, estimated] = bst_vgarch_likelihood(model, x, M, R, S, e, n, toFull, toHalf, false);
% for m = 1:nParameters
%   ground(m) = ...
%     bst_vgarch_likelihood(model, x + ((1:nParameters)==m)' * 1e-6, M, R, S, e, n, toFull, toHalf, false) ...
%     - ...
%     bst_vgarch_likelihood(model, x - ((1:nParameters)==m)' * 1e-6, M, R, S, e, n, toFull, toHalf, false);
% end
% figure(1); clf(1); set(1, 'WindowStyle', 'docked'); plot(1:nParameters, estimated, 'rx', 1:nParameters, ground/2e-6, 'bo');
% 
% ground = zeros(nParameters); [~, ~, estimated] = bst_vgarch_likelihood(model, x, M, R, S, e, n, toFull, toHalf, false);
% for r = 1:nParameters
%   for c = 1:nParameters
%     ground(r,c) = ...
%       bst_vgarch_likelihood(model, x + ((1:nParameters)==r)' * 1e-6 + ((1:nParameters)==c)' * 1e-6, M, R, S, e, n, toFull, toHalf, false) ...
%       - bst_vgarch_likelihood(model, x + ((1:nParameters)==r)' * 1e-6 - ((1:nParameters)==c)' * 1e-6, M, R, S, e, n, toFull, toHalf, false) ...
%       - bst_vgarch_likelihood(model, x - ((1:nParameters)==r)' * 1e-6 + ((1:nParameters)==c)' * 1e-6, M, R, S, e, n, toFull, toHalf, false) ...
%       + bst_vgarch_likelihood(model, x - ((1:nParameters)==r)' * 1e-6 - ((1:nParameters)==c)' * 1e-6, M, R, S, e, n, toFull, toHalf, false);
%   end
% end
% figure(2); clf(2); set(2, 'WindowStyle', 'docked'); plot(1:nParameters*nParameters, estimated(:), 'rx', 1:nParameters*nParameters, ground(:)/(2e-6*2e-6), 'bo');

%% Optimization
if method < 3
  %% quadratic programming & preconditioned conjugate gradient descent
  
  % Default options
  options = optimset('Diagnostics', 'off', 'GradObj', qualifier, ...
    'TolFun', tolerance, 'TolX', tolerance, 'TolCon', tolerance, ...
    'MaxFunEvals', 1e12, 'MaxIter', maxIterations);
  
  % Choice of descent method
  if method == 0 % quadratic programming
    options.Algorithm = 'active-set';
  elseif method == 1 % preconditioned conjugate gradient descent
    options.Algorithm = 'interior-point';
  elseif method == 2 % semiquadratic programming
    options.Algorithm = 'sqp';
  end
  
  % Flag to display results
  if VERBOSE
    options.Display = 'iter';
  else
    options.Display = 'off';
  end
  
  % Initialize the tracking vectors
  if nOutputs > 12
    stepsizes = zeros(maxIterations+1, 1);
    iterationTimes = zeros(maxIterations+1,1);
    logLikelihoods = zeros(maxIterations+1, 1);
    options.OutputFcn = @bst_vgarch_fmincon_iteration;
  end
  
  % fmincon with active-set uses quadratic programming
  if model
    nParameters = M*(M+1)/2 + R*M*M + S*M*M;
  else
    nParameters = M*(M+1)/2 + R*M*(M+1)/2*M*(M+1)/2 + S*M*(M+1)/2*M*(M+1)/2;
  end
  
  [current, logLikelihood, ~, output] = fmincon( ...
    @(syed) bst_vgarch_likelihood(model, syed, M, R, S, e, n, toFull, toHalf, true), ... % optimization function (true to do sign change for minimization)
    [W(:); C(:); D(:)], ... % initial values
    [], [], ... % no linear inequality constraints
    [], [], ... % no linear equality constraints
    [0; -Inf(nParameters-1, 1)], Inf(nParameters, 1), ... % only bound W_11^2 to be positive because of sign ambiguity
    @(syed) deal(bst_vgarch_stability(model, syed, M, R, S, toFull, 0.99998), []), ... % stability: sum of all C's and D's should have no eigenvalue greater than 1
    options);
  
  % final assignment
  if model
    W = current(1:M*(M+1)/2); % convert homoscedastic term into covariance by squaring
    C = reshape(current(M*(M+1)/2 + (1:R*M*M)), M, M, R); % only reshaping is necessary for conditional heteroscedasticity coefficients
    D = reshape(current(M*(M+1)/2 + R*M*M + (1:S*M*M)), M, M, S); % only reshaping is necessary for autoregression in variance coefficients
  else % for the vec model, the coefficients just need to be reshaped into interpretable parts
    W = current(1:M*(M+1)/2); % homoscedastic term, unvectorized
    C = reshape(current(M*(M+1)/2 + (1:R*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, R*M*(M+1)/2); % conditional heteroscedasticity term, unvectorized
    D = reshape(current(M*(M+1)/2 + R*M*(M+1)/2*M*(M+1)/2 + (1:S*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, S*M*(M+1)/2); % autoregression in variance, unvectorized
  end
  if nOutputs > 5 % to also get the information matrix, we have to call bst_vgarch_likelihood again
    [logLikelihood, ~, information] = bst_vgarch_likelihood(model, [W(:); C(:); D(:)], M, R, S, e, n, toFull, toHalf, false);
  else % avoid an extra call to the likelihood function by using the result from fmincon
    logLikelihood = -1 * logLikelihood;    
  end
  iterations = output.iterations;
  
  % for PCG, the iteration count is off by 1 because it doesn't count the initial step
  if method == 1
    iterations = iterations + 1;
  end
  
  % there is no "primal" or "dual" gap so set them to empty
  rho = -1; primalResidual = []; dualResidual = [];
  
elseif method == 3
  %% alternating-direction method of multipliers - full covariance matrix substitution
  
  % store the past to help the computation later
  past = reshape(e(:, fliplr(hankel(n(end) - (1:R), fliplr(n-R)))), R*M*(M+1)/2, N);
  
  rho = ceil(rho); % make an integer for easy display
  
  % Initialize substitutions & Lagrange multipliers
  H = repmat(reshape(W(toFull),M,M), [1 1 N]); % substituted covariance matrices
  eta = zeros(M,M,N); % Lagrange multilpiers corresponding to the covariance matrices
  invEye = repmat(eye(M), [1 1 N]); % used to take the inverse of the covariance matrices in direction 2
  
  % Keep initial values for use when updating rho
  initialW = W; initialC = C;
  
  % Initialize coefficients to be used in direction 1 of the covariance substitution approaches:
  %    _                                               _ 
  %   | 1 e_{(R+1)-1}^T e_{(R+1)-2}^T ... e_{(R+1)-R}^T |
  %   | 1 e_{(R+2)-1}^T e_{(R+2)-2}^T ... e_{(R+2)-R}^T |
  %   | 1 e_{(R+3)-1}^T e_{(R+3)-2}^T ... e_{(R+3)-R}^T |
  %   | .       .             .       .         .       | \otimes P (where P is the half-vectorization of a matrix that is 1 off- and sqrt(1/2) on-diagonal
  %   | .       .             .        .        .       |
  %   | .       .             .         .       .       |
  %   | 1 e_{(N-1)-1}^T e_{(N-1)-2}^T ... e_{(N-1)-R}^T |
  %   |_1 e_{    N-1}^T e_{    N-2}^T ... e_{    N-R}^T_|
  coefficients = [ones(N, 1) zeros(N, R*M*(M+1)/2)]; % coefficient matrix for estimating parameters given covariances
  for r = 1:R % for each lag in heteroscedasticity, add the corresponding lag column in the coefficients
    coefficients(:,1 + (r-1)*M*(M+1)/2 + (1:M*(M+1)/2)) = e(:, n-r)';
  end % by preallocation above, the 1s are already in the matrix
  [U,lambda,V] = svd(coefficients, 'econ'); % economic SVD decomposition to save space & memory
  P = ones(M) + (sqrt(1/2) - 1) * eye(M); P = P(tril(true(M))); % weights for weighted least-squares: 1 off diagonal, 1/sqrt(2) on diagonal
  coefficients = kron(V * diag(1 ./ diag(lambda)) * U', diag(1./P)); % pseudo-inverse of coefficients saved onto itself for ease
  P = repmat(P, N, 1); % weights to multiply to vech(H - N) in direction 1
  
  % Initialize primal & dual residuals for all iterations (for iteration 0, assume the "old" values used in the dual residual are zero)
  primalResidual = zeros(maxIterations+1, 1);
  dualResidual = zeros(maxIterations+1, 1);
  primalResidual(1) = sqrt( sum(sum(sum( (reshape(bst_bsxfun(@plus, W(toFull), C(toFull,:)*past),M,M,N) - H).^2 ))) ) / ( sqrt(N) * M*M );
  dualResidual(1) = sqrt(rho^2 * ...
    sum(( H(toHalf) ).^2) ...
    + sum(sum(sum( bst_bsxfun(@times, reshape(H(toHalf), [M*(M+1)/2 1 N]), reshape(e(:,n), [1 M*(M+1)/2 N])).^2 ))) ...
  ) / ( N * (M*M + M*M*M*M) );

  % Initialize iteration specifics if desired
  fixedIterations = zeros(maxIterations+1, 1);
  if nOutputs > 12
    stepsizes = zeros(maxIterations+1, 1);
    iterationTimes = zeros(maxIterations+1,1);
    logLikelihoods = zeros(maxIterations+1, 1);
  end

  % adjust tolerance if using a final Newton step
  if strcmp(qualifier, 'on')
    tolerance = tolerance*10;
  end
  
  % Initialize display (and zeroth iteration)
  if VERBOSE
    % set up structure and use the iteraiton function I have for fmincon (improves consistency)
    starting = bst_vgarch_likelihood(model, [W; C(:); D(:)], M, R, S, e, n, toFull, toHalf, false);
    optimValues = struct( ...
      'iteration', 0, ...
      'fval', -1 * starting, ...
      'stepsize', sum([W(:); C(:); D(:)].^2) ...
    );
    bst_vgarch_fmincon_iteration([W; C(:); D(:)], optimValues, 'init'); % use init part of the iteration function since we havent started the loop yet
        
    % current status
    fprintf(1, '%5s %19s %19s %19s %19s %10s %s\n', 'count', 'loglikelihood', 'change', 'primal gap', 'dual gap', 'rho', 'parameters');
    fprintf(1, '%5i %19.12f %19s %19s %19s %10i %5s %5s %s\n', 0, -1 * optimValues.fval, '---', '---', '---', rho, '---', '---', sprintf(' %12.9f', [W(:); C(:); D(:)]));

    % plot initial primal & dual gaps
    if nOutputs > 17
      dLikelihood = zeros(maxIterations+1, 1);
      dLikelihood(1) = optimValues.fval;
      gapHandle = figure; clf(gapHandle); set(gapHandle, 'windowstyle', 'docked'); ax = axes('Parent', gapHandle);
      
      [ax, hLike, hPrimal] = plotyy(ax, ...
        0:maxIterations, log10(tolerance)*ones(1, maxIterations+1), ...
        0:maxIterations, log10(tolerance)*ones(1, maxIterations+1) ...
      );
      delete(hLike);
      set(ax, 'fontname', 'Times New Roman', 'fontsize', 12);
      set(hPrimal, 'LineWidth', 2, 'Color', [0 0 0], 'LineStyle', ':', 'Marker', 'none');
      
      set(ax(1), 'xlim', [-1 1], 'ylim', [log10(tolerance)-2 1], 'ytick', log10(tolerance):0);
      hold(ax(1), 'on');
      hLike = plot(ax(1), 0:maxIterations, log10(dLikelihood), 'b');
      hLikeBad = plot(ax(1), find(dLikelihood < 0)-1, log10(abs(dLikelihood(dLikelihood < 0))), 'bx', 'MarkerSize', 16, 'LineWidth', 2);
      hold(ax(1), 'off');
      
      set(ax(2), 'xlim', [-1 1], 'ylim', [log10(tolerance)-2 1], 'ytick', log10(tolerance):0);
      hold(ax(2), 'on');
      hPrimal = plot(ax(2), 0:maxIterations, log10(primalResidual), 'r');
      hDual = plot(ax(2), 0:maxIterations, log10(dualResidual), 'Color', [0 0.5 0]);
      hold(ax(2), 'off');
      
      set(hPrimal, 'LineWidth', 4, 'Marker', 'o'); set(hDual, 'LineWidth', 3, 'Marker', 'o'); set(hLike, 'LineWidth', 1, 'Marker', 'o');
      legend(ax(1), [hPrimal hDual hLike], 'primal gap', 'dual gap', '\Delta likelihood', 'Location', 'Southwest');
      xlabel(ax(1), 'iteration'); ylabel(ax(1), 'log-\Delta likelihood'); ylabel(ax(2), 'gap size'); title(ax(1), ['\rho = ' int2str(rho)]);
      print(1, '-dpng', '-r144', ['~/Desktop/garch/gap_covarianceFull_' int2str(M) '_rho' int2str(rho) '_' sprintf('%0.3i', 0) '.png']);
    end

    % save the old values for use in the next iteration
    oldW = W; oldC = C; oldD = D;
  end
  
  % Iterate until epsilon convergence
  counter = 1; % counter for rho that is too low
  for iterations = 1:maxIterations-1

    % Direction 2: H given w & C (through \tilde{H})
    tildeH = reshape(bst_bsxfun(@plus, W(toFull), C(toFull,:) * past), M, M, N) + eta; % \tilde{H} for direction 2 & primal gap
    oldH = H; % save the old values for use in the gap calculation
    bad = true(1,N); % indices that require gradient descent after the fixed-point iterations

%     tic
%     gradientDescent = zeros(M,M,N); likelihood = cell(N,1); likelihoodFull = cell(N,1);
%     for idx = 1:N
%       best = oldH(:,:,idx);
%       likelihood{idx} = @(syed) log(det(reshape(syed(toFull),M,M))) ...
%           + trace(reshape(syed(toFull),M,M)\eMatrix(:,:,n(idx))) ...
%           + rho * sum(sum(( reshape(syed(toFull),M,M) - tildeH(:,:,idx) ).^2));
%       likelihoodFull{idx} = @(syed) log(det(syed)) + trace(syed\eMatrix(:,:,n(idx))) + rho * sum(sum(( syed - tildeH(:,:,idx) ).^2));
%       best = fminunc( ...
%         likelihood{idx}, ...
%         best(tril(true(M))), ...
%         optimset('display', 'off', 'LargeScale', 'off'));
%       gradientDescent(:,:,idx) = reshape(best(toFull),M,M);
%     end
%     toc

    % Direction 2: fixed-point iterations using MMX
    updateH = H; % to hold the fixed-point updates
    previous = zeros(1,50); % tracker = [N zeros(1,5000)];
    while fixedIterations(iterations+1) < 50 && sum(bad) % first 100
      invH = mmx('backslash', H(:,:,bad), invEye(:,:,bad));
      updateH(:,:,bad) = tildeH(:,:,bad) + 1/(2 * rho) * mmx('mult', invH, mmx('mult', eMatrix(:,:,n(bad)) - H(:,:,bad), invH)); % fixed-point iteration
      bad(bad) = max(max(abs( H(:,:,bad) - updateH(:,:,bad) ))) > tolerance*M/rho; % first test to see where there is any change
      if all(~bad)
        break
      end
      H(:,:,bad) = updateH(:,:,bad); % update those for which there is some change
      fixedIterations(iterations+1) = fixedIterations(iterations+1)+1;
      previous(fixedIterations(iterations+1)) = sum(bad);
      % tracker(fixedIterations(iterations+1)+1) = sum(bad);
    end
    while sum(bad) && sum(bad) < previous(1) % until only bad guys are left
      invH = mmx('backslash', H(:,:,bad), invEye(:,:,bad));
      updateH(:,:,bad) = tildeH(:,:,bad) + 1/(2 * rho) * mmx('mult', invH, mmx('mult', eMatrix(:,:,n(bad)) - H(:,:,bad), invH)); % fixed-point iteration
      bad(bad) = max(max(abs( H(:,:,bad) - updateH(:,:,bad) ))) > tolerance*M/rho; % first test to see where there is any change
      if all(~bad)
        break
      end
      H(:,:,bad) = updateH(:,:,bad); % update those for which there is some change
      fixedIterations(iterations+1) = fixedIterations(iterations+1)+1;
      previous = [previous(2:end) sum(bad)];
      % tracker(fixedIterations(iterations+1)+1) = sum(bad);
    end
    % tracker = tracker(1:fixedIterations(iterations+1)+1);

    % Direction 2: fixed-point iteration cleanup via gradient descent
    % -- If you want to avoid all the fixed-point stuff, comment out the fixed-point iteration loop above only (it then works normally)
    if any(bad) && sum(bad) < 100
      for idx = find(bad)
        best = oldH(:,:,idx);
        try
          best = fminunc( ...
            @(syed) 2*log(det(reshape(syed(toFull),M,M))) ...
              + trace(reshape(syed(toFull),M,M)*eMatrix(:,:,n(idx))) ...
              + rho * sum(sum(( reshape(syed(toFull),M,M) - tildeH(:,:,idx) ).^2)), ...
            best(tril(true(M))), ...
            optimset('display', 'off', 'LargeScale', 'off'));
          H(:,:,idx) = reshape(best(toFull),M,M);
        catch ME %#ok<NASGU> if gradient descent fails, just use the old value and hope the next iteration fixes everything
          H(:,:,idx) = best;
        end
      end
    end

    % primal & dual gaps
    primalGap = tildeH - eta - H; % gap for covariance matrix constraints: \tilde{H} = hcev^{-1} (w + C e_{n-1}) + N and we want hcev^{-1} (w + C e_{n-1}) - H
    primalResidual(iterations+1) = sqrt( sum( primalGap(toHalf).^2 ) ) / ( N * sqrt(M*(M+1)/2) );
    dualResidual(iterations+1) = sqrt(rho^2 * ...
      sum(( H(toHalf)-oldH(toHalf) ).^2) ...
      + sum(sum(sum( bst_bsxfun(@times, reshape(H(toHalf)-oldH(toHalf), [M*(M+1)/2 1 N]), reshape(e(:,n), [1 M*(M+1)/2 N])).^2 ))) ...
    ) / ( N * sqrt(M*(M+1)/2 + M*(M+1)/2*M*(M+1)/2) );
    % first term is change in h, second term is change in h weighted by the signal (so (\hat{h}_n - h_n)^T * e_{n-1} to be exact)

    % First (and only) break: a gap has closed. I use 4x tolerance because the final iterations are slow and have little benefit.
    if primalResidual(iterations+1) < tolerance && dualResidual(iterations+1) < tolerance
      break;
    elseif counter > 2 && ( ...
      (primalResidual(iterations+1) > tolerance && primalResidual(iterations+1) > primalResidual(iterations)*1.05) ...
      || ...
      (dualResidual(iterations+1) > tolerance && dualResidual(iterations+1) > dualResidual(iterations)*1.05) ...
    )
      if counter > maxIterations/10 % If after 10% iterations you are diverging, you're actually pretty close so just stop the algorithm
        break;
      end
      rho = ceil(rho*1.1); % increase rho
      W = initialW; C = initialC; eta = zeros(M,M,N); H = repmat(reshape(W(toFull),M,M), [1 1 N]); % reset the initializations
      counter = 1; % reset the counter so that the new rho gets 3 iterations
      continue
    end
    counter = counter+1; % counter for number of iterations for this choice of rho
    % It may be useful to do a one-step gradient search to find the best value after this approach, which would get around the 3x tolerance issue

    % Dual update - I should think of a way to check that this is right ...
    eta = eta + primalGap; % for the multipliers corresponding to the precision matrix constraints

    % Direction 1: w, C given G
    % this is the only time where the model makes a difference: BEKK requires a regression where C is a Kronecker product of a matrix with itself
    % weighted linear regression: P * w + \sum_{r=1}^R P{ C_r e_{n-r} = h_n - eta_n, P is diagonal matrix with 1s except at (i,i) points where it is 1/sqrt(2)
    initial = coefficients * (P .* (H(toHalf) - eta(toHalf))); % use the pre-computed pseudo-ivnerse to get the weighted linear estimation of w & C
    W = initial(1:M*(M+1)/2); % homoscedastic covariance (only the lower half-triangle)
    C = reshape(initial(M*(M+1)/2 + (1:R*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, R*M*(M+1)/2); % the conditional heteroscedasticity coefficients (C * vech(e))
    
    % Display
    if VERBOSE
      oldLikelihood = -1 * optimValues.fval; % old likelihood value is computed in previous iteration
      
      % set up structure and use the iteraiton function I have for fmincon (improves consistency)
      optimValues = struct( ...
        'iteration', iterations, ...
        'fval', bst_vgarch_likelihood(model, [W; C(:); D(:)], M, R, S, e, n, toFull, toHalf, true), ...
        'stepsize', sum([W(:) - oldW(:); C(:) - oldC(:); D(:) - oldD(:)].^2) ...
        );
      bst_vgarch_fmincon_iteration([W(:); C(:); D(:)], optimValues, 'iter'); % use iter part of the iteration function since we're in the loop
      
      % Current status
      if mod(iterations, 30) == 0
        disp(' ');
        fprintf(1, '%5s %19s %19s %19s %19s %10s %s\n', 'count', 'loglikelihood', 'change', 'primal gap', 'dual gap', 'rho', 'parameters');
      end
      fprintf(1, '%5i %19.12f %19.12f %19.10f %19.10f %10i %5i %5i %s\n', ...
        iterations, ...
        -1 * optimValues.fval, ...
        -1 * optimValues.fval - oldLikelihood, ...
        primalResidual(iterations+1), ...
        dualResidual(iterations+1), ...
        rho, ...
        fixedIterations(iterations+1), sum(bad), ...
        sprintf(' %12.9f', [W(:); C(:); D(:)]));
      
      % plot current progress in terms of change in log-likelihood as well as current primal & dual gaps
      if nOutputs > 17
        dLikelihood(iterations+1) = -1 * optimValues.fval - oldLikelihood;
        delete(hPrimal); delete(hDual); delete(hLike); delete(hLikeBad);
        
        hold(ax(1), 'on');
        hLike = plot(ax(1), 0:iterations, log10(abs(dLikelihood(1:iterations+1))), 'b');
        hLikeBad = plot(ax(1), find(dLikelihood < 0)-1, log10(abs(dLikelihood(dLikelihood < 0))), 'bx', 'MarkerSize', 16, 'LineWidth', 2);
        hold(ax(1), 'off');
        
        hold(ax(2), 'on');
        hPrimal = plot(ax(2), 0:iterations, log10(primalResidual(1:iterations+1)), 'r');
        hDual = plot(ax(2), 0:iterations, log10(dualResidual(1:iterations+1)), 'Color', [0 0.5 0]);
        hold(ax(2), 'off');
        
        set(ax, 'xlim', [-1 iterations+1]);
        set(hPrimal, 'LineWidth', 4, 'Marker', 'none'); set(hDual, 'LineWidth', 3, 'Marker', 'none'); set(hLike, 'LineWidth', 1, 'Marker', 'none');
        
        if iterations < 100
          set([hPrimal hDual hLike], 'Marker', 'o');
        elseif iterations == 100
          legend(ax(1), [hPrimal hDual hLike], 'primal gap', 'dual gap', '\Delta likelihood', 'Location', 'Southwest');
        end
        if iterations > 50 && all(primalResidual(iterations+1 + (-20:0)) < 1e-4) && all(dualResidual(iterations+1 + (-20:0)) < 1e-4)
          set(ax(1), 'ylim', [-13 -4], 'ytick', -12:-5);
          set(ax(2), 'ylim', [log10(tolerance)-3 -4], 'ytick', log10(tolerance):0.2:-4.2);
        end
        
        drawnow
        print(1, '-dpng', '-r144', ['~/Desktop/garch/gap_covarianceFull_' int2str(M) '_rho' int2str(rho) '_' sprintf('%0.3i', iterations) '.png']);
        % move the images to a big server like dynamics.usc.edu and then run
        % convert -delay 20 ~/Desktop/garch/gap_2_rho10_*.png -delay 60 $(ls -1 ~/Desktop/garch/gap_2_rho10_*.png | tail -n 1) -loop 0 ~/Desktop/garch/gap_2_rho10.gif
      end

      % save the old values for use in the next iteration
      oldW = W; oldC = C; % oldD = D;
    end
    
  end

  % final iteration: if desired, we perform two Newton steps with a stepsize found by setting a parabola to the loglikelihood in the search direction
  if strcmp(qualifier, 'on')
    for final = 1:M-1
      % final iteration: get the direction to perform our one line search
      [logLikelihood, score, information] = bst_vgarch_likelihood(model, [W(:); C(:); D(:)], M, R, S, e, n, toFull, toHalf, false); % particulars at final point
      direction = (information \ score); % preconditioned gradient descent direction

      % final iteration: the marker "behind" the current point (based on the direction) for interpolation
      stepBack = -8; % get a wide interval for the parabolic interpolation
      back = bst_vgarch_likelihood(model, [W(:); C(:); D(:)] + stepBack*direction, M, R, S, e, n, toFull, toHalf, false); % value for marker "behind" final point
      while (isinf(back) || isnan(back)) && stepBack < -1e-6 % while the marker "behind" is infeasible but still far away
        stepBack = stepBack/2; % step the marker closer towards the final point 
        back = bst_vgarch_likelihood(model, [W(:); C(:); D(:)] + stepBack*direction, M, R, S, e, n, toFull, toHalf, false); % value for marker "behind" final point
      end

      % final iteration: the marker "in front" of the current point (based on the direction) for interpolation
      stepFront = 8; % get a wide interval for the parabolic interpolation
      front = bst_vgarch_likelihood(model, [W(:); C(:); D(:)] + stepFront*direction, M, R, S, e, n, toFull, toHalf, false); % value for marker "in front" of final point
      while (isinf(front) || isnan(front)) && stepFront > 1e-6 % while the marker "in front" is infeasible but still far away
        stepFront = stepFront/2; % step the marker closer towards the final point
        front = bst_vgarch_likelihood(model, [W(:); C(:); D(:)] + stepFront*direction, M, R, S, e, n, toFull, toHalf, false); % value for marker "in front" of final point
      end

      % final iteration: if admissible, use the two markers to build a parabola and use its peak as the stepsize
      if stepBack < -1e-6 && stepFront > 1e-6 % both markers are far enough away that interpolation is a good idea

        % parabolic fit by linear estimation of the coefficients
        % estimated = [stepBack^2 stepBack 1; 0 0 1; stepFront^2 stepFront 1] \ [back; logLikelihood; front]; % three points makes a parabola; least-squares to fit
        % best = min(max(-estimated(2)/(2*estimated(1)), stepBack), stepFront); % maximum will be at -b/2a for ax^2 + bx + c, but we also control for feasibility

        % parabolic fit by direct computation of the quadratic and linear coefficients
        % quadratic = (stepBack*(front-logLikelihood) - stepFront*(back-logLikelihood)) / (stepBack * stepFront * (stepFront-stepBack));
        % linear = (stepFront^2*(back-logLikelihood) - stepBack^2*(front-logLikelihood)) / (stepBack * stepFront * (stepFront-stepBack));
        % best = min(max(-linear/(2*quadratic), stepBack), stepFront); % maximum will be at -b/2a for ax^2 + bx + c, but we also control for feasibility

        % parabolic fit AND computation of the best stepsize (so that I don't have to divide by the normalization factor
        best = min(max( ...
          (stepBack^2*(front-logLikelihood) - stepFront^2*(back-logLikelihood)) ...
          / 2 / ... maximum will be at -b/2a for ax^2 + bx + c
          (stepBack*(front-logLikelihood) - stepFront*(back-logLikelihood)) ...
          , stepBack), stepFront); % we also control for stability by projecting back onto the interval

        % code to test the parabolic fit; it's pretty damn good
        % steps = -3:0.1:3;
        % loglikelihoodCurve = zeros(length(steps),1);
        % for iStep = 1:length(steps)
        %   loglikelihoodCurve(iStep) = bst_vgarch_likelihood(model, [W(:); C(:); D(:)] + steps(iStep) * (information \ score), M, R, S, e, n, toFull, toHalf, false);
        % end
        % figure(1); clf(1); set(1, 'WindowStyle', 'docked');
        % bestLL = bst_vgarch_likelihood(model, [W(:); C(:); D(:)] + best * direction, M, R, S, e, n, toFull, toHalf, false);
        % plot(steps, loglikelihoodCurve, 'kx', steps, estimated(1)*steps.^2 + estimated(2)*steps + estimated(3), 'r-', best, bestLL, 'bo', 'MarkerSize', 12, 'LineWidth', 2);

        W = W + best * direction(1:M*(M+1)/2); % move the homoscedastic term based on the parabolic interpolation in the preconditioned direction
        C = C + best * reshape(direction(M*(M+1)/2 + (1:R*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, R*M*(M+1)/2); % also move conditonal heteroscedasticity coefficietns
        D = D + best * reshape(direction(M*(M+1)/2 + R*M*(M+1)/2*M*(M+1)/2 + (1:S*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, S*M*(M+1)/2); % and autoregression-in-variance coefficients

        % final iteration: update the residuals being tracked
        iterations = iterations+1;
        primalResidual(iterations+1) = sqrt( sum( primalGap(toHalf).^2 ) ) / sqrt( N * M*(M+1)/2 );
        dualResidual(iterations+1) = sqrt(rho^2 * ...
          sum(( H(toHalf)-oldH(toHalf) ).^2) ...
          + sum(sum(sum( bst_bsxfun(@times, reshape(H(toHalf)-oldH(toHalf), [M*(M+1)/2 1 N]), reshape(e(:,n), [1 M*(M+1)/2 N])).^2 ))) ...
        ) / ( N * (M*(M+1)/2 + M*(M+1)/2*M*(M+1)/2) );
      end

    end
    
  end
  
  % final iteration: update the other trackers being used
  if VERBOSE
    oldLikelihood = -1 * optimValues.fval; % old likelihood value is computed in previous iteration

    % set up structure and use the iteraiton function I have for fmincon (improves consistency)
    optimValues = struct( ...
      'iteration', iterations, ...
      'fval', bst_vgarch_likelihood(model, [W; C(:); D(:)], M, R, S, e, n, toFull, toHalf, true), ...
      'stepsize', sum([W(:) - oldW(:); C(:) - oldC(:); D(:) - oldD(:)].^2) ...
      );
    bst_vgarch_fmincon_iteration([W(:); C(:); D(:)], optimValues, 'iter'); % use iter part of the iteration function since we're in the loop

    % Current status
    if mod(iterations, 30) == 0
      disp(' ');
      fprintf(1, '%5s %19s %19s %19s %19s %10s %s\n', 'count', 'loglikelihood', 'change', 'primal gap', 'dual gap', 'rho', 'parameters');
    end
    fprintf(1, '%5i %19.12f %19.12f %19.10f %19.10f %10i %5.0d %5.0d %s\n', ...
      iterations, ...
      -1 * optimValues.fval, ...
      -1 * optimValues.fval - oldLikelihood, ...
      primalResidual(iterations+1), ...
      dualResidual(iterations+1), ...
      rho, ...
      stepBack, stepFront, ...
      sprintf(' %12.9f', [W(:); C(:); D(:)]));

    % plot current progress in terms of change in log-likelihood as well as current primal & dual gaps
    if nOutputs > 17
      dLikelihood(iterations+1) = -1 * optimValues.fval - oldLikelihood;
      delete(hPrimal); delete(hDual); delete(hLike); delete(hLikeBad);

      hold(ax(1), 'on');
      hLike = plot(ax(1), 0:iterations, log10(abs(dLikelihood(1:iterations+1))), 'b');
      hLikeBad = plot(ax(1), find(dLikelihood < 0)-1, log10(abs(dLikelihood(dLikelihood < 0))), 'bx', 'MarkerSize', 16, 'LineWidth', 2);
      hold(ax(1), 'off');

      hold(ax(2), 'on');
      hPrimal = plot(ax(2), 0:iterations, log10(primalResidual(1:iterations+1)), 'r');
      hDual = plot(ax(2), 0:iterations, log10(dualResidual(1:iterations+1)), 'Color', [0 0.5 0]);
      hold(ax(2), 'off');

      set(ax, 'xlim', [-1 iterations+1]);
      set(hPrimal, 'LineWidth', 4, 'Marker', 'none'); set(hDual, 'LineWidth', 3, 'Marker', 'none'); set(hLike, 'LineWidth', 1, 'Marker', 'none');

      if iterations < 100
        set([hPrimal hDual hLike], 'Marker', 'o');
      elseif iterations == 100
        legend(ax(1), [hPrimal hDual hLike], 'primal gap', 'dual gap', '\Delta likelihood', 'Location', 'Southwest');
      end
      if iterations > 50 && all(primalResidual(iterations+1 + (-20:0)) < 1e-4) && all(dualResidual(iterations+1 + (-20:0)) < 1e-4)
        set(ax(1), 'ylim', [-13 -4], 'ytick', -12:-5);
        set(ax(2), 'ylim', [log10(tolerance)-3 -4], 'ytick', log10(tolerance):0.2:-4.2);
      end

      drawnow
      print(1, '-dpng', '-r144', ['~/Desktop/garch/gap_covarianceFull_' int2str(M) '_rho' int2str(rho) '_' sprintf('%0.3i', iterations) '.png']);
      % move the images to a big server like dynamics.usc.edu and then run
      % convert -delay 20 ~/Desktop/garch/gap_2_rho10_*.png -delay 60 $(ls -1 ~/Desktop/garch/gap_2_rho10_*.png | tail -n 1) -loop 0 ~/Desktop/garch/gap_2_rho10.gif
    end

    % save the old values for use in the next iteration
    oldW = W; oldC = C; % oldD = D;
  end  

  % Display final result
  if VERBOSE
    fprintf(1, '%5s %19.12f %19.12f %19.10f %19.10f %10i %s\n', ...
      '---->', ...
      logLikelihood, ...
      logLikelihood-starting, ...
      primalResidual(iterations+1), ...
      dualResidual(iterations+1), ...
      rho, ...
      sprintf(' %12.9f', [W(:); C(:); D(:)]));
  elseif nOutputs > 5
    [logLikelihood, ~, information] = bst_vgarch_likelihood(model, [W(:); C(:); D(:)], M, R, S, e, n, toFull, toHalf, false); % particulars at final point
  end
  
  % trim primal & dual residual vectors
  if nOutputs > 9
    primalResidual = primalResidual(1:iterations+1);
    dualResidual = dualResidual(1:iterations+1);
  end
  
elseif method == 4
  %% Bregman-proximal trust-region method (Chretien2012) -- NOT STARTED
  
else % ERRUR
  error('BrainStorm:bst_vgarch:method', 'method should be 0, 1, 2 or 3.');
end

%% Cleanup: runtime & clip the tracking vectors (stepsizes, iteration times, log-likelihoods, and possibly handles for the moving GIFs)
if nOutputs > 12
  stepsizes = stepsizes(1:iterations+2);
  iterationTimes = iterationTimes(1:iterations);
  logLikelihoods = logLikelihoods(1:iterations);
  if exist('handles', 'var')
    handles = handles(1:iterations);
  end
end

% convert to the friendlier full-vectorization model
% W = reshape(W(toFull), M, M); % convert W to a symmetric covariance matrix
% if isempty(model) || ~model
%   
%   C = reshape(C, M*(M+1)/2, M*(M+1)/2, R); C = C(toFull, toFull, :); % fill in C with the coefficients caused by the symmetry
%   C(:, reshape(logical(1-eye(M)), M*M, 1), :) = C(:, reshape(logical(1-eye(M)), M*M, 1), :)/2; % coefficient distribution to the full model
%   % the coefficients corresponding to the off-diagonal elements of e are all divided by 2, since they were smashed together by symmetry in the estimation
%   
%   D = reshape(D, M*(M+1)/2, M*(M+1)/2, S); D = D(toFull, toFull, :); % fill in D with the coefficients caused by the symmetry
%   D(:, reshape(logical(1-eye(M)), M*M, 1), :) = D(:, reshape(logical(1-eye(M)), M*M, 1), :)/2; % coefficient distribution to the full model
%   % the coefficients corresponding to the off-diagonal elements of e are all divided by 2, since they were smashed together by symmetry in the estimation
%   
% end

% runtime
elapsed = toc(totalStart);
if VERBOSE
  fprintf(1, 'Elapsed time is %0.6f seconds.\n', elapsed);
end

%% Nested function: iteration at each fmincon if needed
function stop = bst_vgarch_fmincon_iteration(x, optimValues, state)
  
  W2 = sum(x(1:M*(M+1)/2).^2);
  if model
    C2 = 2 * sum(x(M*(M+1)/2 + (1:R*M*M)).^2)^2;
    % D2 = 2 * sum(x(M*(M+1)/2 + (1:S*M*M)).^2)^2;
  else
    C2 = sum(x(M*(M+1)/2 + (1:R*M*(M+1)/2*M*(M+1)/2)).^2);
    % D2 = sum(x(M*(M+1)/2 + (1:S*M*(M+1)/2*M*(M+1)/2)).^2);
  end
  
  if strcmp(state, 'init') % Starting graph: initial point
    
    if ~isempty(ax1)
      hold(ax1, 'on');
      handles(1) = plot(ax1, C2, W2); % plot initial point
      set(handles(1), 'Color', lineColor, 'LineWidth', 2, 'Marker', 'o', 'MarkerSize', 8, 'MarkerFaceColor', lineColor, 'MarkerEdgeColor', lineColor);
      hold(ax1, 'off'); title(ax1, qualifier);      
      oldW2 = W2; oldC2 = C2; % oldD2 = D2; % save this iteration for the next loop
    end
    
  elseif strcmp(state, 'iter') % Next graph: iteration
    
    % progress for this iteration
    if nOutputs > 12
      if optimValues.iteration > 0
        stepsizes(optimValues.iteration+1) = optimValues.stepsize; % Displacement (same as stepsize for other algorithms)
        % for ADMM, the stepsize is calculated in the construction of optimValues before a call to this function
      end
      iterationTimes(optimValues.iteration+1) = toc(totalStart); % Iteration time
      if nOutputs > 14
        logLikelihoods(optimValues.iteration+1) = -1 * optimValues.fval; % Current loglikelihood
        % for ADMM, the likelihood is calculated in the construction of optimValues before a call to this function
      end
    end
    
    if ~isempty(ax1) % Path tracing
      hold(ax1, 'on');
      handles(optimValues.iteration+1) = plot(ax1, [oldC2 C2], [oldW2 W2]); % plot path
      set(handles(optimValues.iteration+1), 'Color', lineColor, 'LineWidth', 2); % make sure it's thick
      delete(handles(1)); % remove the previous big circle
      handles(1) = plot(ax1, C2, W2); % plot current location with a new big circle
      set(handles(1), 'Color', lineColor, 'LineWidth', 2, 'Marker', 'o', 'MarkerSize', 8, 'MarkerFaceColor', lineColor, 'MarkerEdgeColor', lineColor); % make sure it's thick and visible
      hold(ax1, 'off'); title(ax1, [qualifier ' descent']); % add a title
      drawnow; % I want it shown before it is saved as a figure
      saveas(hLikelihood, [handlePrefix '_' sprintf('%0.4i', length(regexp(strtrim(ls('-1', [handlePrefix '_*.fig'])), '\n', 'split'))) '.fig'], 'fig'); % save as figure because I can load them later to ensure I like the fonts, etc.
      oldW2 = W2; oldC2 = C2; % oldD2 = D2; % save this iteration for the next loop
    end
    
  end
  
  stop = false; % Do not stop the optimization; let fmincon stop it for us if it wants to
  % this value is useless for ADMM, so ADMM calls this function without looking at the return value
  
end

end

%% ===================================================================== LIKELIHOOD MODEL ======================================================================
function [logLikelihood, score, information, h] = bst_vgarch_likelihood(model, coefficients, M, R, S, e, n, toFull, toHalf, flagMinimize)
% BST_VGARCH_LIKELIHOOD   Calculation of multivariate GARCH log-likelihood
%                         (and gradients if asked for)
%
% Inputs:
%   model               - model used to compute vector GARCH log-likelihood
%   coefficients        - coefficients corresponding to the model above
%                         for model = 'BEKK', we have coefficients as the
%                         one-dimensional stack of the vectorized matrices
%                           W: the square root of the homoscedastic covariance
%                              (so WW^T is the homoscedastic covariance)
%                           C: the conditional heteroscedasticity coefficients
%                              (which look like C e_{n-1} e_{n-1}^T C^T)
%                           D: the autoregression-in-variance coefficients
%                              (which look like D h_{n-1} h_{n-1}^T D^T)
%                         for model = 'vec', we have coefficients as the
%                         one-dimensional stack of the vectorized matrices
%                           W: the homoscedastic covariance
%                           C: the conditional heteroscedasticity coefficients
%                           D: the autoregression-in-variance coefficients
%   M                   - dimensionality of signals (innovations)
%   R                   - number of lags in the conditional heteroscedasticity
%   S                   - number of lags in the autoregression in variance
%   e                   - Outer product of signal, vectorized
%   n                   - indices for "present" values
%                         (based on number of trials and maximum lag)
%   toFull              - indexing from the lower half-triangle vectorization
%                         to the original symmetric matrix
%                         --> used to reconstruct W
%   toHalf              - indexing from the original symmetric matrix to the
%                         lower half-triangle vectorization
%   flagMinimize        - flag to invert values if minimizing
%                         (e.g. if using for fmincon or other standard package)
% 
% Outputs:
%   loglikelihood       - Average log-likelihood over number of samples
%   score               - Average score over number of samples
%   information         - Average information matrix over number of samples
%                         (so expected value of Hessian)

%% Setup
N = length(n);
W = coefficients(1:M*(M+1)/2); % homoscedastic term in half-vectorized form
if model % assign and then convert to vec model
  
  % assignment of BEKK parameters
  CBEKK = reshape(coefficients(M*(M+1)/2 + (1:R*M*M)), M, M, R); % turn into matrices for the Kronecker mapping to a vec model
  DBEKK = reshape(coefficients(M*(M+1)/2 + R*M*M + (1:S*M*M)), M, M, S); % turn into matrices for the Kronecker mapping to a vec model
  
  % conversion to vec model
  C = zeros(M*M, M*M, R); D = zeros(M*M, M*M, S);
  for r = 1:R
    C(:, (r-1)*M*M + (1:M*M)) = kron(CBEKK(:,:,r), CBEKK(:,:,r)); % for the conditionally heteroscedastic term, the vectorization is a kronecker product
  end
  for s = 1:S
    D(:, (s-1)*M*M + (1:M*M)) = kron(DBEKK(:,:,r), DBEKK(:,:,r)); % as is true for the autoregression-in-variance term
  end
  
  % conversion to half-vectorized model
  C = reshape(C(toHalf, toHalf, :), M*(M+1)/2, R*M*(M+1)/2); % keep coefficients corresponding to lower half-triangle of e & h, and then stack in block row
  D = reshape(D(toHalf, toHalf, :), M*(M+1)/2, S*M*(M+1)/2); % keep coefficients corresponding to lower half-triangle of e & h, and then stack in block row
  
else
  C = reshape(coefficients(M*(M+1)/2 + (1:R*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, R*M*(M+1)/2); % conditional heteroscedasticity term, unvectorized
  D = reshape(coefficients(M*(M+1)/2 + R*M*(M+1)/2*M*(M+1)/2 + (1:S*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, S*M*(M+1)/2); % autoregression in variance, unvectorized
end
flagGradient = (nargout > 1); % check to see if we also want to send back the gradients
if flagGradient % weight vector to compute the gradient from the half-vectorized model
  gradientWeights = ones(M) - 1/2 * eye(M);
  gradientWeights = gradientWeights(tril(true(M)));
  score = zeros(M*(M+1)/2+R*M*(M+1)/2*M*(M+1)/2+S*M*(M+1)/2*M*(M+1)/2, 1);
  information = zeros(M*(M+1)/2+R*M*(M+1)/2*M*(M+1)/2+S*M*(M+1)/2*M*(M+1)/2);
  eHalf = e; % need the half-vectorized innovations to compute the gradient & information matrix
end

%% Variances

% ARCH portion of model
h = zeros(M*(M+1)/2, size(e,2));
h(:, setdiff(1:size(e,2), n)) = repmat(W, 1, size(e,2)-N); % set homoscedastic term to be the default variance for the first few timepoints in each trial

% vectorize the computation of the covariance half-vectors
h(:,n) = bst_bsxfun(@plus, W, C * reshape(e(:, hankel(n(1)-(1:R), n-R)), R*M*(M+1)/2, N)); % use bsxfun because w is a column addition to each timepoint column
if ~isempty(D) % generalized addendum requires a for loop because it is feedback
  for idx = n
    h(:,idx) = h(:,idx) + D * reshape(h(:, idx-(1:S)), S*M*(M+1)/2, 1);
  end
end

% un-vectorized: unnecessary because the vectorized solution is stable
% for idx = n
%   h(:,idx) = W + C * reshape(e(:, idx-(1:R)), R*M*(M+1)/2, 1) + D * reshape(h(:, idx-(1:S)), S*M*(M+1)/2, 1);
% end

%% Outputs: log-likelihood, score & information matrix

% this part requires the linear algebra of the full symmetric forms, so we have to convert the variances & estimates to that form
H = reshape(h(toFull, :), M, M, size(e,2)); % convert to full symmetric matrix form for inversion
e = reshape(e(toFull, :), M, M, size(e,2)); % convert to full symmetric vectorized form for marginal likelihood estimation

% vectorize the computation of the log-likelihood -- very UNSTABLE
% e = reshape(e, [M M size(e,2)]); % convert to full symmetric vectorized form for marginal likelihood estimation
% determinant = mmx('chol', H(:,:,n), []);
% determinant = prod(reshape(determinant(repmat(logical(eye(M)), [1 1 N])), M, N)).^2;
% exponential = mmx('backslash', H(:,:,n), e(:,:,n));
% exponential = sum(reshape(exponential(repmat(logical(eye(M)), [1 1 N])), M, N));
% logLikelihood = -N * M*log(2*pi) - sum(log(determinant)) - sum(exponential);
% if flagGradient % DO NOT USE - the gradient has not been mapped to the half-vectorized coefficients in the vec model, or modified for the BEKK model
%   invH = mmx('backslash', H(:,:,n), repmat(eye(M), [1 1 N]));
%   projection = mmx('mult', invH, mmx('mult', H(:,:,n) - e(:,:,n), invH));
%   score = [ ...
%     reshape(projection, M*M, N); ...
%     reshape(reshape(projection, [M*M 1 N]) .* reshape(e(:,:,n-1), 1, M*M, N), M*M*M*M, N); ...
%     zeros(M*M*M*M*S, 1) ...
%   ];
% 
%   % =========== HERE: map the score as a fcn of hcev(w) and hcev(C \times \cdot) to a fcn of w and C ===========
% 
%   % aggregate marginal gradients to get information matrix & score
%   information = score * score';
%   score = -sum(score, 2);
% end

% un-vectorized computation of the log-likelihood, which is required to get a solution for high M
logLikelihood = -N * M*log(2*pi); % constant term (we put the 1/2 in at the end)
for idx = n % sum over all timepoints

  % log-likelihood
  try % when the matrix is positive definite, the fastest determinant calculation is a Cholesky decomposition 
    logLikelihood = logLikelihood - 2*log(prod(diag(chol(H(:,:,idx))))) - sum(diag(H(:,:,idx) \ e(:,:,idx)));
  catch %#ok<CTCH> % when the matrix is not positive definite, the likelihood is ZERO and thus the log-likelihood is minus infinity
    logLikelihood = -Inf;
    if flagGradient % also we set the gradient & information matrix to infinity to account for this egregious behavior by calling functions
      score = Inf;
      information = Inf;
    end
    break; % no need to compute the rest of the terms, so go straight to the break
  end
  % - log(|H_n|) - e_n^T H_n^{-1} e_n = - log(|H_n|) - Tr{H_n^{-1} e_n e_n^T}
  % the determinant is the square of the product of the diagonal entries in the Cholesky decomposition

  % do not calculate the gradient unless it is asked of the calling program (i.e. when an analytic gradient is desired)
  if flagGradient

    % marginal gradient
    projection = ( H(:,:,idx) \ (e(:,:,idx) - H(:,:,idx)) ) / H(:,:,idx); % the key term: H_n^{-1} ( H_n - e_n e_n^T ) H_n
    marginal = [ ...
      gradientWeights .* projection(tril(true(M))); ... % gradient with respect to w for time n is just the projection term times -1/2
      reshape((gradientWeights .* projection(tril(true(M)))) * reshape(eHalf(:, idx-(1:R)), [], 1)', [], 1); ... % gradient with respect to C for time n is vectorized form of projection * [e_{n-1}; e_{n-2}; ...; e_{n-r}] transposed
      zeros(S*M*(M+1)/2*M*(M+1)/2, 1) ... % for now, ignore autoregression in variance because it is iterative and annoying
    ];

    % aggregation to form score (and information matrix): because the signal is Markov, the gradient is a sum of these (marginal) gradients over timepoints
    score = score + marginal; % the marginal gradient is added
    information = information + marginal * marginal'; %  % the marginal information is added

  end

end

%% Scale by number of timepoints to get rates (and I think increase robustness)
logLikelihood = logLikelihood / (N*2) / M; % divide by 2 to account for the 1/2 term, divide by N and M to get a better dynamic range value for the optimization
if flagGradient
  score = score / N / M; % divide by N to get the expectation, and divide by M to get a better dynamic range value for the optimization
  information = information / N; % divide by N to get the expectation
  % note that we do NOT divide the information matrix by M^2, since we need the raw information for causality testing (via the Wald statistic)
end

%% Flip sign of log-likelihood and score so minimization is optimal point
if flagMinimize 
  logLikelihood = -logLikelihood;
  if flagGradient
    score = -score;
  end
end

end

%% ============================================================== INITIALIZATION OF COEFFICIENTS ===============================================================
function [W, C, D] = bst_vgarch_initialize(model, e, n, M, R, S, N, toFull, toHalf, method, qualifier, maxIterations, tolerance, VERBOSE)
% BST_VGARCH_INITIALIZE   Initial estimate of vector GARCH model fit to the
%                         signal with chosen order in conditional
%                         heteroscedasticity and autoregression-in-variance.
% Inputs:
%   model         - Model used for vector GARCH modelling
%                   'BEKK', 'bekk': use symmetric, constrained BEKK model
%                    'vec',  'vec': use more general, simpler vec model
%                   ^-- by default, the 'vec' model is used
%                   [Z: string]
%   e             - Vector signal that is autoregressive in variance
%                   [E: M x NT double]
%   n             - Samples for the "present" timepoints
%                   [T: (N-max(R,S))T x 1 positive integer]
%   R             - Order of autoregression-in-variance
%                   [R: scalar nonnegative integer]
%   S             - Order of heteroscedasticity-in-variance
%                   [S: scalar nonnegative integer]
%                   ^-- currently, only S = 0 works
%   method        - Method used to generate initialization
%                   0: variance targeting w/ active-set method (same as SQP)
%                   1: variance targeting w/ interior point method (PCG)
%                   2: variance targeting w/ SQP (quadratic programming)
%                   3: least-squares using single-sample variance estimation
%                   X: error
%                   [C: positive integer]
%   qualifier     - Desired modified for the method chosen above in C
%                   For PCG & SQP,
%                   'off': use numerical gradient
%                    'on': use analytic gradient
%                   For least-squares estimation, the qualifier is not used.
%                   [Q: string]
%   maxIterations - Maximum number of iterations for any of the methods
%                   [I: scalar positive integer]
% Outputs:
%   W             - Homoscedastic covariance term in GARCH model
%                   [W: M(M+1)/2 x 1 positive double]
%   C             - Conditional heteroscedasticity coefficients in GARCH model
%                   [C: M x M x R or M(M+1)/2 x M(M+1)/2 x R positive double]
%   D             - Autoregression-in-variance coefficients in GARCH model
%                   [D: M x M x S or M(M+1)/2 x M(M+1)/2 x S positive double]

% @=============================================================================
% This software is part of the Brainstorm software:
% http://neuroimage.usc.edu/brainstorm
% 
% Copyright (c)2000-2013 Brainstorm by the University of Southern California
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html.
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Syed Ashrafulla, 2012

% This code uses mmx, a software package that allows many matrix multiplications
% at the same time. To do this, you need to download build_mmx from
% http://www.mathworks.com/matlabcentral/fileexchange/37515-mmx-multithreaded-matrix-operations-on-n-d-matrices
% and then go to that folder and edit build_mmx.m to modify one line and add two lines:
% -- on line 80, change the iteration to go from 2 to 3
% ---- do not do this if you have the Intel MKL installed. You won't have this
% ---- installed unless you know it immediately because you are a developer.
% -- on line 196, add "syed = pwd"
% -- on line 214, add "cd(syed)"
% Next, you need to update /Applications/MATLABR2013a.app/bin/mexopts.sh. Use
% vim and the command
% :%s/10.7/10.8
% and you will get 5 substitutions. Write it out using :wq! to force overwriting
% and then go back into MatLab and hit mex -setup. Pick 1 to get your edits onto
% your use of MatLab.
% Finally, build_mmx will work with 1 stupid type conversion warning!

if method < 3
  %% Scalar variance targeting
  
  % Default options
  options = optimset('Diagnostics', 'off', 'GradObj', qualifier, ...
    'TolFun', tolerance, 'TolX', tolerance, 'TolCon', tolerance, ...
    'MaxFunEvals', 1e12, 'MaxIter', maxIterations);
  
  % Choice of descent method
  if method == 0 % quadratic programming
    options.Algorithm = 'active-set';
  elseif method == 1 % preconditioned conjugate gradient descent
    options.Algorithm = 'interior-point';
  elseif method == 2 % semiquadratic programming
    options.Algorithm = 'sqp';
  end
  
  % Flag to display results
  if VERBOSE
    options.Display = 'iter';
  else
    options.Display = 'off';
  end
  
  % from MFEToolbox
%   startingOptions = optimset('fminunc');
%   startingOptions.LargeScale = 'off';
%   startingOptions.Display = 'none';
%   [initial,~,~,intercept] = scalar_vt_vech(reshape(e(toFull,n), M, M, N),[],R,0,S,[],[],startingOptions);
%   testW = intercept(tril(true(M))); % the homoscedastic covariance's lower half-triangle
%   if model
%     testC = kron(sqrt(initial(1:R)), eye(M)); % CMFE_r = sqrt(initial(r)) * I
%     testD = kron(sqrt(initial(R+(1:S))), eye(M)); % DMFE_s = sqrt(initial(R+s)) * I
%   else
%     testC = kron(initial(1:R), eye(M*(M+1)/2)); % in vec form, C_r = CMFE_r^T \otimes CMFE_r from vectorization; since CMFE_r = sqrt(initial(r)) * I, C_r is initial(r) * I
%     testD = kron(initial(R+(1:S)), eye(M*(M+1)/2)); % in vec form, D_s = DMFE_r^T \otimes DMFE_r from vectorization; since DMFE_s = sqrt(initial(R+s)) * I, D_s is initial(R+s) * I
%   end

  % initialize by guessing
  covariance = sum(e(:,n),2)/(N-(n(1)-1)); % normalize by number of timepoints minus number of trials for unbiased covariance
  best = bst_vgarch_initialize_likelihood(zeros(R+S,1), M, R, S, e, n, toFull, toHalf, covariance, true); % first grid point: 0 conditional heteroscedasticity
  C = zeros(R,1); D = zeros(S,1); % will hold the best coefficients after the rough grid search
  testC = ones(R,1)*[0 .005 .01 .025 .05]; % possible choices for initialization (start with a small conditional heteroscedasticity parameter)
  if S == 0; testD = zeros(0,1); else testD = ones(S,1)*[.99 .95 .90 .5 0]; end; % if no autoregression-in-variance, ignore D; otherwise, test high values for D
  for iC = 1:size(testC,2)
    for iD = 1:size(testD,2)
      current = bst_vgarch_initialize_likelihood([testC(:,iC); testD(:,iD)], M, R, S, e, n, toFull, toHalf, covariance, true); % next grid point
      if current < best % save if better
        best = current; C = testC(:,iC); D = testD(:,iD);
      end
    end
  end
  
  % fmincon with active-set uses quadratic programming
  current = fmincon( ...
    @(syed) bst_vgarch_initialize_likelihood(syed, M, R, S, e, n, toFull, toHalf, covariance, true), ... % optimization function (sign changed)
    [C; D], ... % initial values
    [], [], ... % no linear inequality constraints
    [], [], ... % no linear equality constraints
    [], [], ... % only bound W_11^2 to be positive because of sign ambiguity
    @(syed) deal(sum(syed) - 0.99998, []), ... % stability: sum of all C's and D's should not be bigger than 1
    options);

  % assignment
  W = (1 - sum(current))*sum(e(:,n),2)/N; % the homoscedastic covariance's lower half-triangle
  if model
    C = kron(sqrt(current(1:R).'), eye(M)); % CMFE_r = sqrt(initial(r)) * I
    if S > 0
      D = kron(sqrt(current(R+(1:S)).'), eye(M)); % DMFE_s = sqrt(initial(R+s)) * I
    else
      D = zeros(M, 0);
    end
  else
    C = kron(current(1:R).', eye(M*(M+1)/2)); % in vec form, C_r = CMFE_r^T \otimes CMFE_r from vectorization; since CMFE_r = sqrt(initial(r)) * I, C_r is initial(r) * I
    if S > 0
      D = kron(current(R+(1:S)).', eye(M*(M+1)/2)); % in vec form, D_s = DMFE_r^T \otimes DMFE_r from vectorization; since DMFE_s = sqrt(initial(R+s)) * I, D_s is initial(R+s) * I
    else
      D = zeros(M*(M+1)/2, 0);
    end
  end

  % use this to show that MFEToolbox & us are the same
%   disp([testW' initial'; W' current'; [W' current'] - [testW' initial']]);
  
elseif model
  %% Least-squares on the first estimate of each diagonal entry (so a diagonal BEKK model initialization)
  
  W = zeros(M); C = zeros(M, M, R); D = zeros(M, M, S);
  for m = 1:M
    
    % regression assuming h_n^2 = e_n^2
    regressor = [ones(1,N); zeros(R, N)];
    for r = 1:R
      regressor(r+1, :) = e((m-1)*M+m,n-r);
    end
    initial = e((m-1)*M+m,n) / regressor;
    
    % assignment: homoscedastic for first element, heteroscedastic for the remaining elements, and nearly-zero for autoregression-in-variance
    W(m,m) = initial(1); C(m, m, :) = initial(1+(1:R)); D(m, m, :) = initial((1+R+(1:S)));
  end
  
  W = W(tril(true(M))); % only save the lower half-triangle of W
  C = reshape(C, M, R*M); D = reshape(D, M, S*M); % turn the other two variables into a block matrix
  
else
  %% Least-squares on the first estimate of the covariance matrices (which would be the vec formulation)

  % linear regression: I * w + \sum_{r=1}^R C_r e_{n-r} = e_n
  % the linear form is [I (e_{n-1}^T \otimes I) (e_{n-2}^T \otimes I) ... (e_{n-R}^T \otimes I)] [w; vec(C_1); vec(C_2); ...; vec(C_R)] = e_n
  %  _ stack time over columns, so that it is                                            _   _            _     _       _
  % | I (e_{(R+1)-1}^T \otimes I) (e_{(R+1)-2}^T \otimes I) ... (e_{(R+1)-R}^T \otimes I) | |      w       | = | e_{R+1} |
  % | I (e_{(R+2)-1}^T \otimes I) (e_{(R+2)-2}^T \otimes I) ... (e_{(R+2)-R}^T \otimes I) | |   vec(C_1)   | = | e_{R+2} |
  % | I (e_{(R+3)-1}^T \otimes I) (e_{(R+3)-2}^T \otimes I) ... (e_{(R+3)-R}^T \otimes I) | |   vec(C_2)   | = | e_{R+3} |
  % | .             .                         .             .               .           ) | |      .       | = |    .    |
  % | .             .                         .              .              .           ) | |      .       | = |    .    |
  % | .             .                         .               .             .           ) | |      .       | = |    .    |
  % | I (e_{(N-1)-1}^T \otimes I) (e_{(N-1)-2}^T \otimes I) ... (e_{(N-1)-R}^T \otimes I) | | vec(C_{R-1}) | = | e_{N-1} |
  % |_I (e_{    N-1}^T \otimes I) (e_{    N-2}^T \otimes I) ... (e_{    N-R}^T \otimes I)_| |_  vec(C_R)  _| = |_e_{  N}_|
  % the final note is that the matrix on the left can be turned into one Kronecker product
  %   | 1 e_{(R+1)-1}^T e_{(R+1)-2}^T ... e_{(R+1)-R}^T |             |      w       | = | e_{R+1} |
  %   | 1 e_{(R+2)-1}^T e_{(R+2)-2}^T ... e_{(R+2)-R}^T |             |   vec(C_1)   | = | e_{R+2} |
  %   | 1 e_{(R+3)-1}^T e_{(R+3)-2}^T ... e_{(R+3)-R}^T |             |   vec(C_2)   | = | e_{R+3} |
  %   | .       .             .       .         .       |             |      .       | = |    .    |
  % ( | .       .             .        .        .       | \otimes I ) |      .       | = |    .    |
  %   | .       .             .         .       .       |             |      .       | = |    .    |
  %   | 1 e_{(N-1)-1}^T e_{(N-1)-2}^T ... e_{(N-1)-R}^T |             | vec(C_{R-1}) | = | e_{N-1} |
  %   |_1 e_{    N-1}^T e_{    N-2}^T ... e_{    N-R}^T_|             |_  vec(C_R)  _| = |_e_{  N}_|
  coefficients = [ones(N, 1) zeros(N, R*M*(M+1)/2)]; % build the coefficient matrix: first, the column of ones
  for r = 1:R  % for each lag in heteroscedasticity, add the corresponding lag column in the coefficients
    coefficients(:,1 + (r-1)*M*(M+1)/2 + (1:M*(M+1)/2)) = e(:, n-r)';
  end
  warning('off', 'MATLAB:rankDeficientMatrix');
  initial = ...
    kron(coefficients, eye(M*(M+1)/2)) ... % kronecker product described above (corresponding to vectorization of C in solution)
    \ reshape(e(:, n), [], 1); % then hit the least-squares estimation to get the resulting vector of initial starting points
  warning('on', 'MATLAB:rankDeficientMatrix');
  
  % assignment
  W = initial(1:M*(M+1)/2); % homoscedastic covariance (only the lower half-triangle)
  C = reshape(initial(M*(M+1)/2 + (1:R*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, R*M*(M+1)/2); % the conditional heteroscedasticity coefficients (C * vech(e))
  D = repmat(eye(M*(M+1)/2), [1 S]) * 0.01; % the autoregression-in-variance coefficients (D * vech(e))
  
end

end

%% ========================================================== LIKELIHOOD MODEL FOR VARIANCE TARGETING ==========================================================
function [logLikelihood, score, information, h] = bst_vgarch_initialize_likelihood(coefficients, M, R, S, e, n, toFull, toHalf, covariance, flagMinimize)
% BST_VGARCH_INITALIZELIKELIHOOD  Calculation of scalar multivariate GARCH
%                                 log-likelihood, or equivalently the vector
%                                 GARCH log-likelihood when every coefficient
%                                 matrix is a constant times an identity matrix
%                                 (including gradients if desired).
%
% Inputs:
%   model               - model used to compute vector GARCH log-likelihood
%   coefficients        - coefficients corresponding to the model above
%                         for model = 'BEKK', we have coefficients as the
%                         one-dimensional stack of the vectorized matrices
%                           W: the square root of the homoscedastic covariance
%                              (so WW^T is the homoscedastic covariance)
%                           C: the conditional heteroscedasticity coefficients
%                              (which look like C e_{n-1} e_{n-1}^T C^T)
%                           D: the autoregression-in-variance coefficients
%                              (which look like D h_{n-1} h_{n-1}^T D^T)
%                         for model = 'vec', we have coefficients as the
%                         one-dimensional stack of the vectorized matrices
%                           W: the homoscedastic covariance
%                           C: the conditional heteroscedasticity coefficients
%                           D: the autoregression-in-variance coefficients
%   M                   - dimensionality of signals (innovations)
%   R                   - number of lags in the conditional heteroscedasticity
%   S                   - number of lags in the autoregression in variance
%   e                   - Outer product of signal, vectorized
%   n                   - indices for "present" values
%                         (based on number of trials and maximum lag)
%   toFull              - indexing from the lower half-triangle vectorization
%                         to the original symmetric matrix
%                         --> used to reconstruct W
%   toHalf              - indexing from the original symmetric matrix to the
%                         lower half-triangle vectorization
%   flagMinimize        - flag to invert values if minimizing
%                         (e.g. if using for fmincon or other standard package)
% 
% Outputs:
%   loglikelihood       - Average log-likelihood over number of samples
%   score               - Average score over number of samples
%   information         - Average information matrix over number of samples
%                         (so expected value of Hessian)

%% Setup
N = length(n);
W = (1 - sum(coefficients))*covariance; % homoscedastic term in half-vectorized form, computed from covariance identitiy: R_ee = C + a * R_ee + b * R_ee
C = kron(coefficients(1:R).', eye(M*(M+1)/2)); % conditionally heteroscedastic term
D = kron(coefficients(R+(1:S)).', eye(M*(M+1)/2)); % autoregression-in-variance term

flagGradient = (nargout > 1); % check to see if we also want to send back the gradients
if flagGradient 
  gradientWeights = ones(M) - 1/2 * eye(M); gradientWeights = gradientWeights(tril(true(M))); % weight vector to compute the gradient from the half-vectorized model
  score = zeros(M*(M+1)/2+R*M*(M+1)/2*M*(M+1)/2+S*M*(M+1)/2*M*(M+1)/2, 1);
  information = zeros(M*(M+1)/2+R*M*(M+1)/2*M*(M+1)/2+S*M*(M+1)/2*M*(M+1)/2);
  eHalf = e; % need the half-vectorized innovations to compute the gradient & information matrix
end

%% Variances

% ARCH portion of model
h = zeros(M*(M+1)/2, size(e,2));
h(:, setdiff(1:size(e,2), n)) = repmat(W, 1, size(e,2)-N); % set homoscedastic term to be the default variance for the first few timepoints in each trial

% vectorize the computation of the covariance half-vectors
h(:,n) = bst_bsxfun(@plus, W, C * reshape(e(:, hankel(n(1)-(1:R), n-R)), R*M*(M+1)/2, N)); % use bsxfun because w is a column addition to each timepoint column
if ~isempty(D) % generalized addendum requires a for loop because it is feedback
  for idx = n
    h(:,idx) = h(:,idx) + D * reshape(h(:, idx-(1:S)), S*M*(M+1)/2, 1);
  end
end

%% Outputs: log-likelihood, score & information matrix

% this part requires the linear algebra of the full symmetric forms, so we have to convert the variances & estimates to that form
H = reshape(h(toFull, :), M, M, size(e,2)); % convert to full symmetric matrix form for inversion
e = reshape(e(toFull, :), M, M, size(e,2)); % convert to full symmetric vectorized form for marginal likelihood estimation

% vectorize the computation of the log-likelihood -- very UNSTABLE
determinant = mmx('chol', H(:,:,n), []);
determinant = prod(reshape(determinant(repmat(logical(eye(M)), [1 1 N])), M, N)).^2;
exponential = mmx('backslash', H(:,:,n), e(:,:,n));
exponential = sum(reshape(exponential(repmat(logical(eye(M)), [1 1 N])), M, N));
logLikelihood = -N * M*log(2*pi) - sum(log(determinant)) - sum(exponential);
if flagGradient % DO NOT USE - the gradient has not been mapped to the half-vectorized coefficients in the vec model, or modified for the BEKK model
  invH = mmx('backslash', H(:,:,n), repmat(eye(M), [1 1 N]));
  projection = mmx('mult', invH, mmx('mult', H(:,:,n) - e(:,:,n), invH));
  score = [ ...
    bst_bsxfun(@times, gradientWeights, reshape(projection(toHalf), M*(M+1)/2, N)); ...
    squeeze(sum(bst_bsxfun(@times, reshape(projection(toHalf), [M*(M+1)/2 1 N]), reshape(eHalf(:, hankel(n(1)-(1:R), n-R)), M*(M+1)/2, R, N)))); ...
    zeros(S, N) ...
  ];

  % aggregate marginal gradients to get information matrix & score
  information = score * score';
  score = -sum(score, 2);
end

%% Scale by number of timepoints to get rates (and I think increase robustness)
logLikelihood = logLikelihood / (N*2) / M; % divide by 2 to account for the 1/2 term, divide by N and M to get a better dynamic range value for the optimization
if flagGradient
  score = score / N / M; % divide by N to get the expectation, and divide by M to get a better dynamic range value for the optimization
  information = information / N; % divide by N to get the expectation
  % note that we do NOT divide the information matrix by M^2, since we need the raw information for causality testing (via the Wald statistic)
end

%% Flip sign of log-likelihood and score so minimization is optimal point
if flagMinimize 
  logLikelihood = -logLikelihood;
  if flagGradient
    score = -score;
  end
end

end

%% ================================================================= CONSTRAINTS (BEKK & VEC) ==================================================================
function largest = bst_vgarch_stability(model, coefficients, M, R, S, toFull, worst)
% BST_VGARCH_STABILITY    Test stability of current coefficients depending on
%                         the type of multivariate GARCH model
%
% Inputs:
%   model               - model used to compute vector GARCH log-likelihood
%   coefficients        - coefficients corresponding to the model above
%                         for model = 'BEKK', we have coefficients as the
%                         one-dimensional stack of the vectorized matrices
%                           W: the square root of the homoscedastic covariance
%                              (so WW^T is the homoscedastic covariance)
%                           C: the conditional heteroscedasticity coefficients
%                              (which look like C e_{n-1} e_{n-1}^T C^T)
%                           D: the autoregression-in-variance coefficients
%                              (which look like D h_{n-1} h_{n-1}^T D^T)
%                         for model = 'vec', we have coefficients as the
%                         one-dimensional stack of the vectorized matrices
%                           W: the homoscedastic covariance
%                           C: the conditional heteroscedasticity coefficients
%                           D: the autoregression-in-variance coefficients
%   M                   - dimension of signal
%   R                   - number of lags in the conditional heteroscedasticity
%   S                   - number of lags in the autoregression in variance
%   worst               - largest value allowed in eigenvalue
%                         (choose a value near but not 1, to enforce stability)
% Outputs:
%   largest             - distance from max eigenvalue to largest allowed
%                         -- will be negative whenever an eigenvalue exceeds
%                         -- will be positive when all eigenvalues are below

%% Coefficient reformatting for stability calculation

% W does not change the stability criterion
% W = current(1:M*(M+1)/2); W = W(toFull); % homoscedastic term, converted to the full symmetric matrix (but still vectorized)

if model % assign and then convert to vec model
  
  % assignment of BEKK parameters
  CBEKK = reshape(coefficients(M*(M+1)/2 + (1:R*M*M)), M, M, R); % turn into matrices for the Kronecker mapping to a vec model
  DBEKK = reshape(coefficients(M*(M+1)/2 + R*M*M + (1:S*M*M)), M, M, S); % turn into matrices for the Kronecker mapping to a vec model
  
  % conversion to vec model
  C = zeros(M*M, M*M, R); D = zeros(M*M, M*M, S);
  for r = 1:R
    C(:, :, r) = kron(CBEKK(:,:,r), CBEKK(:,:,r)); % for the conditionally heteroscedastic term, the vectorization is a kronecker product
  end
  for s = 1:S
    D(:, :, s) = kron(DBEKK(:,:,r), DBEKK(:,:,r)); % as is true for the autoregression-in-variance term
  end
  
else
  C = reshape(coefficients(M*(M+1)/2 + (1:R*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, M*(M+1)/2, R); % conditional heteroscedasticity term, unvectorized
  D = reshape(coefficients(M*(M+1)/2 + R*M*(M+1)/2*M*(M+1)/2 + (1:S*M*(M+1)/2*M*(M+1)/2)), M*(M+1)/2, M*(M+1)/2, S); % autoregression in variance, unvectorized
  
  % conversion to the full-vectorized model
  C = C(toFull, toFull, :); D = D(toFull, toFull, :);
end

%% Stability metric
% stability test: are all the eigenvalues within the unit circle? use worst to prevent quasistability

largest = abs(eig(sum(C, 3) + sum(D, 3))) - worst;
% largest must be < 0, meaning all eigenvalues of sum(C_r) + sum(D_s) are less than worst (usually some value close to but not quite 1)

end
