R/FitMultivariateGarch.R
In R-Finance/Meucci: Collection of functionality ported from the MATLAB code of Attilio Meucci.
#' Estimation of multivariate GARCH models
#'
#' @param returns : [matrix] (T x N) returns so rows must correspond to time and columns to assets
#' @param demean : [scalar] specifies whether returns should be demeaned (if demean = 1) or not to estimate the model; default value is 1.
#' @param eps : [scalar] used in enforcing a_ii + b_ii <= 1 - eps; the default value is zero
#' @param df : [scalar] degree of freedom for the t-distribution; the default value is 500 to make it, basically, normal
#'
#' @return mu : [vector]
#' @return ATMF : [matrix] coefficient matrix A-tilde (in the notation of the paper)
#' @return BTMF : [matrix] coefficient matrix B-tilde (in the notation of the paper)
#' @return CTMF : [matrix] coefficient matrix C-tilde (in the notation of the paper)
#' @return Hhat : [matrix] forecasted conditional covariance matrix
#'
#' @note Code for MATLAB initially written by Olivier Ledoit and Michael Wolf
#'
#' @references
#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170},
#' "E 136 - Equity market: multivariate GARCH process".
#'
#' See Meucci's script for "FitMultivariateGarch.m"
#'
#' @author Xavier Valls \email{flamejat@@gmail.com}
#' @export
FitMultivariateGarch = function( returns, demean = 1, eps = 0, df = 500 )
{
if (eps < 0)
{
error("eps must be a (small) positive number")
}
# Initialization
T = nrow( returns );
N = ncol( returns );
if ( 1 == demean )
{
mu = matrix(apply( returns, 2, mean ));
returns = returns-mu[ array(1,T ), ];
}
S = t(returns) %*% returns / ( T - 1 );
x = t(returns);
A = matrix( 0, N, N );
B = matrix( 0, N, N );
C = matrix( 0, N, N );
# Rescale Data
scale = sqrt( matrix( apply( t(x)^2, 2, mean ) ) );
x = x / scale[ , array( 1, T ) ];
# Estimation of On-Diagonal Elements
h = matrix( 0, N, T );
for( i in 1 : N )
{
# Likelihood Maximization
q = garch1f4( t( x[ i, ] ), eps, df)$q;
A[ i, i ] = q[ 2 ];
B[ i, i ] = q[ 3 ];
C[ i, i ] = q[ 1 ];
h[ i, ] = filter( cbind( 0, q[ 2] ), cbind( 1, -q[ 3 ] ), x[ i, ]^2 * (df-2) / df, mean( x[ i, ] ^ 2) * ( df-2 ) / df ) +
filter( cbind( 0, q[ 1 ] ), cbind( 1, -q[ 3 ] ), matrix( 1, 1, T) );
}
# First-step Estimation of Off-Diagonal Elements
for( i in 1 : (N-1) )
{
for( j in (i +1) : N )
{
# Likelihood Maximization
theta = garch2f8( x[ i, ] * x[ j, ], C[ i, i ], A[ i, i ], B[ i, i ], x[ i, ]^2, h[ i, ], C[ j, j ], A[ j, j ], B[ j, j ], x[ j, ]^2, h[ j, ], df );
A[ i, j ] = theta[ 2 ];
B[ i, j ] = theta[ 3 ];
C[ i, j ] = theta[ 1 ];
A[ j, i ] = A[ i, j ];
B[ j, i ] = B[ i, j ];
C[ j, i ] = C[ i, j ];
}
}
theta = garch2f8( x[ N, ] * x[ N, ], C[ N, N ], A[ N, N ], B[ N, N ], x[ N, ]^2, h[ N, ], C[ N, N ], A[ N, N ], B[ N, N ], x[ N, ]^2, h[ N, ], df );
A[ N, N ] = theta[ 2 ];
B[ N, N ] = theta[ 3 ];
C[ N, N ] = theta[ 1 ];
A[ N, N ] = A[ N, N ];
B[ N, N ] = B[ N, N ];
C[ N, N ] = C[ N, N ];
# Transformation of Coefficient Matrices
ATMF = minfro( A );
BTMF = minfro( B );
CTMF = minfro( C /( 1 - B )) * ( 1 - BTMF );
# Rescale
#C = C .* (scale * scale');
CTMF = CTMF * ( scale %*% t(scale) );
# Forecast of Conditional Covariance Matrix
Hhat = matrix( 0, N, N );
for( i in 1 : N )
{
for( j in 1 : N )
{
hSeries = filter( cbind( 0, ATMF[ i, j ] ), cbind( 1, -BTMF[ i, j ] ), returns[ , i ] * returns[ , j ], S[ i, j ] ) +
filter( cbind( 0, CTMF[ i, j ] ), cbind( 1, -BTMF[ i, j ] ), matrix( 1, T,1));
Hhat[ i, j ] = hSeries[ T ];
}
}
mu = matrix( mu );
return( list( mu = mu, ATMF = ATMF, BTMF = BTMF, CTMF = CTMF, Hhat = Hhat ) );
}
#' Fit a GARCH(1,1) model with student-t errors
#'
#' @param x : [vector] (T x 1) data generated by a GARCH(1,1) process
#' @param eps : [scalar] used in enforcing a_ii + b_ii <= 1 - eps; the default value is zero
#' @param df : [scalar] degree of freedom for the t-distribution; the default value is 500 to make it, basically, normal
#'
#' @return q : [vector] (4 x 1) parameters of the GARCH(1,1) process
#' @return qerr : [vector] (4 x 1) standard error of parameter estimates
#' @return hf : [scalar] current conditional heteroskedasticity estimate
#' @return hferr : [scalar] standard error on hf
#'
#' @note
#' MATLAB's script initially written by Olivier Ledoit, 4/28/1997
#'
#' Uses a conditional t-distribution with fixed degrees of freedom
#'
#' Difference with garch1f: errors come from the score alone
#'
#' @references
#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170}.
#'
#' See Meucci's script for "FitMultivariateGarch.m"
#'
#' @author Xavier Valls \email{flamejat@@gmail.com}
#' @export
garch1f4 = function( x, eps, df )
{
# Parameters
gold = ( 1 + sqrt(5) ) / 2; # step size increment
tol1 = 1e-7; # for termination criterion
tol2 = 1e-7; # for closeness to boundary
big = 2; # for making the hessian negative definite
maxiter = 50; # maximum number of iterations
n = 30; # number of points on the grid
# Rescale
y = ( array(x) - mean( array(x)))^2;
t = length(y);
scale = sqrt(mean( y^2 ));
y = y / scale;
s = mean(y);
# Grid search
#[ag, bg] = meshgrid( seq( 0, 1 - eps, length = n ), seq( 0, 1-eps, length = n ) );
ag = outer(seq( 0, 1 - eps, length = n )*0, seq( 0, 1-eps, length = n ), FUN = "+" );
bg = outer(seq( 0, 1 - eps, length = n ) , seq( 0, 1-eps, length = n ) * 0, FUN = "+" )
cg = pmax( s * (1-ag-bg), 0 );
likeg = matrix( -Inf, n, n );
for( i in 1 : n )
{
for( j in 1 : (n-i+1) )
{
h = filter( cbind( 0, ag[ i, j ] ), cbind( 1, -bg[ i, j] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, cg[ i, j ]),cbind( 1, - bg[ i, j ] ),matrix(1, t, 1));
likeg[ i, j ] = -sum( log(h) + (df+1) * log( 1 + y /h /df) );
}
}
maxlikeg = max(likeg, na.rm = TRUE);
maxima = which(likeg == maxlikeg); ##ok<MXFND>
maximum = max( maxima );
a = cbind(cg[ maximum ], ag[ maximum], bg[ maximum ]);
best = 0;
da = 0;
#term = 1;
#negdef = 0;
iter = 0;
while( iter < maxiter )
{
iter = iter + 1;
# New parameter1
a = a + gold^best * da;
# Conditional variance
h = filter(cbind( 0, a[ 2 ] ),cbind( 1, -a[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) + filter(cbind( 0, a[ 1 ] ),cbind( 1, -a[ 3 ] ), matrix(1, t, 1 ) );
# Likelihood
if( any( a<0 )||( (a[ 2 ] + a[ 3 ]) > (1 - eps) ))
{
like = -Inf;
} else
{
like = -sum( log( h )+( df + 1 ) * log( 1 + y /h /df) );
}
GG = cbind( filter(cbind(0, 1),cbind( 1, -a[ 3] ), matrix( 1, t, 1 )),
filter(cbind( 0, 1 ),cbind( 1, -a[3] ), y * (df-2) / df ),
filter(cbind( 0, 1 ),cbind( 1, -a[ 3 ]), h ));
colnames(GG)<-NULL;
g1 = ( ( df+1 ) * ( y / ( y + df * h ) ) - 1 ) / h;
G = GG * (g1 %*% matrix( 1, 1, 3));
gra = apply( G, 2, sum );
# Hessian
GG2 = GG[ ,c( 1, 2, 3, 1, 2, 3, 1, 2, 3)] * GG [, c( 1, 1, 1, 2, 2, 2, 3, 3, 3 ) ];
g2 = -((df+1) * ( y /( y + df*h) ) - 1 ) / h^2 - ( df * ( df + 1 ) )*(y / ( y+df * h ) ^2 / h );
HH = matrix( 0, t, 9);
HH[ , 3 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ ,1 ] );
HH[ , 7 ] = HH[ , 3 ];
HH[ , 6 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ , 2 ] );
HH[ , 8 ] = HH[ ,6 ];
HH[ , 9 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ , 3 ] );
H = GG2 * (g2 %*% matrix( 1, 1, 9)) + HH * (g1 %*% matrix( 1, 1, 9)) ;
hes = matrix( apply( H, 2, sum), 3, 3 );
e = eigen(hes);
if( any(e$values > 0) )
{
negdef = 0;
d = min( e$values, max( e$values[ e$values<0 ] ) / big);
hes = e$vectors %*% diag(d, length(d)) %*% t(e$vectors);
} else
{
negdef = 1;
}
# Direction
da = -gra %*% solve( hes );
# Termination criterion
term = (da %*% gra)[1];
if ((term < tol1) && negdef)
{
break;
}
best = 0;
newa = a + gold^(best-1) * da;
if( any(newa < 0 ) || ( newa[ 2 ] + newa[ 3 ] > (1-eps) ) )
{
left = -Inf;
} else
{
h = filter( cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, newa[ 1 ]),cbind( 1, - newa[ 3 ] ),matrix(1, t, 1));
left = -sum( log(h) + (df+1) * log( 1 + y /h /df) )
}
newa = a + gold^best * da;
if( any(newa < 0 ) || ( newa[ 2 ] + newa[ 3 ] > (1-eps) ) )
{
center = -Inf;
} else
{
h = filter( cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, newa[ 1 ]),cbind( 1, - newa[ 3 ] ),matrix(1, t, 1));
center = -sum( log(h) + (df+1) * log( 1 + y /h /df) )
}
newa = a + gold^( best+1 ) * da;
if( any(newa < 0 ) || ( newa[ 2 ] + newa[ 3 ] > (1-eps) ) )
{
right=-Inf;
} else
{
h = filter( cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, newa[ 1 ]),cbind( 1, - newa[ 3 ] ),matrix(1, t, 1));
right = -sum( log(h) + (df+1) * log( 1 + y /h /df) )
}
if( all(like > c( left, center, right)) || all( left > c( center, right )) )
{
while( 1 )
{
best = best-1;
center = left;
newa = a + gold^( best-1 ) * da;
if( any(newa < 0 ) || ( newa[ 2 ] + newa[ 3 ] > (1-eps) ) )
{
left = -Inf;
} else
{
h = filter( cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, newa[ 1 ]),cbind( 1, - newa[ 3 ] ),matrix(1, t, 1));
left = -sum( log(h) + (df+1) * log( 1 + y /h /df) )
}
if( all( center >= c( like, left ) ) )
{
break;
}
}
}else
{
if( all( right > c( left, center ) ) )
{
while( 1 )
{
best = best+1;
center = right;
newa = a + gold^(best+1) *da;
if( any(newa < 0 ) || ( newa[ 2 ] + newa[ 3 ] > (1-eps) ) )
{
right = -Inf;
} else
{
h = filter( cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, newa[ 1 ]),cbind( 1, - newa[ 3 ] ),matrix(1, t, 1));
right = -sum( log(h) + (df+1) * log( 1 + y /h /df) )
}
if( center > right )
{
break;
}
}
}
}
# If stuck at boundary then stop
if( ( center == like ) && ( any(a<tol2) || ( ( a[ 2 ] + a[ 3 ] ) > (1-tol2) ) ) )
{
break;
}
# end of optimization loop
}
if( length(a[ a < tol2 ]) )
{
a[ a < tol2 ] = matrix( 0, 1, length(a[ a<tol2 ]));
}
if( ( a[ 2 ] + a[ 3 ] )> ( 1-tol2 ) )
{
if( a[ 2 ]< ( 1 - tol2 ) )
{
a[ 2 ] = a[ 2 ] + (1 - a[ 2 ]-a[ 3 ] );
}else
{
a[ 3 ] = a[ 3 ] + ( 1 - a[ 2 ]- a[ 3 ] );
}
}
# Estimation error and volatility forecast
#aerr=solve(t(G)%*%G);
tmp = ( t(G) %*% G );
aerr = tmp %*% solve(diag( 1, dim(tmp)));
hf = a[ 1 ] + a[ 2 ] * y[ t ] * (df-2) / df + a[ 3 ] * h[ t ];
gf = c( 1, y[ t ], h[ t ] ) + a[ 3 ] * GG[ t, ];
hferr = gf %*% aerr %*% gf;
aerr = t( diag(aerr) );
# Revert to original scale
a[ 1 ] = a[ 1 ] * scale;
aerr[ 1 ] = aerr[ 1 ] * scale^2;
hf = hf* scale;
hferr = hferr * scale^2;
aerr = sqrt(aerr);
hferr = sqrt(hferr);
q = a;
qerr = aerr;
return( list( q = q, qerr = qerr, hf = hf, hferr = hferr ) );
}
#' Off-diagonal parameter estimation in bivariate GARCH(1,1) when diagonal parameters are given.
#'
#' @param y : [vector] (T x 1) data generated by a GARCH(1,1) process
#' @param c1 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix C
#' @param a1 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix A
#' @param b1 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix B
#' @param y1 : [vector] (T x 1) data generated by a GARCH(1,1) process
#' @param h1 : [vector] (T x 1) data generated by a GARCH(1,1) process
#' @param c2 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix C
#' @param a2 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix A
#' @param b2 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix B
#' @param y2 : [vector] (T x 1) data generated by a GARCH(1,1) process
#' @param h2 : [vector] (T x 1) generated by a GARCH(1,1) process
#' @param df : [scalar] degree of freedom for the t-distribution; the default value is 500 to make it, basically, normal
#'
#' @return q : [vector] (4 x 1) parameters of the GARCH(1,1) process
#' @return qerr : [vector] (4 x 1) standard error of parameter estimates
#' @return hf : [scalar] current conditional heteroskedasticity estimate
#' @return hferr : [scalar] standard error on hf
#'
#' @note
#' MATLAB's code initially written by Olivier Ledoit, 4/28/1997
#'
#' Uses a conditional t-distribution with fixed degrees of freedom
#'
#' Steepest Ascent on boundary, Hessian off boundary, no grid search
#'
#' @references
#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170}.
#'
#' See Meucci's script for "FitMultivariateGarch.m"
#'
#' @author Xavier Valls \email{flamejat@@gmail.com}
#' @export
garch2f8 = function( y, c1, a1, b1, y1, h1, c2, a2, b2, y2, h2, df )
{
# Parameters
gold = ( 1 + sqrt( 5 ) ) / 2; # step size increment
tol1 = 1e-7; # for termination criterion
tol2 = 1e-7; # for closeness to boundary
big = 2; # for making the hessian negative definite
maxiter = 50; # maximum number of iterations
#n=30; # number of points on the grid
# Prepare
t = length( y );
y1 = array( y1 );
y2 = array( y2 );
y = array( y );
s = mean( y );
#s1 = mean(y1);
#s2 = mean(y2);
h1 = array( h1 );
h2 = array( h2 );
# Bounds
low = c(-sqrt(c1*c2), 0, 0 ) + tol2;
high = c( sqrt(c1*c2), sqrt(a1*a2), sqrt(b1*b2) ) - tol2;
# Starting Point
a0 = 0.9 * sqrt( a1 * a2 );
b0 = 0.9 * sqrt( b1 * b2 );
c0 = mean( y ) * ( 1 - a0 - b0 ) %*% ( df - 2 )/df;
c0 = sign( c0 ) * min( abs( c0 ), 0.9 * sqrt( c1*c2 ));
# Initialize optimization
a = cbind( c0, a0, b0 );
best = 0;
da = 0;
#term=1;
#negdef=0;
iter = 0;
while( iter < maxiter )
{
iter = iter + 1;
# New parameter1
a = a + gold^best * da;
# Conditional variance
h = filter(cbind( 0, a[ 2 ] ),cbind( 1, -a[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) + filter(cbind( 0, a[ 1 ] ),cbind( 1, -a[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
# Likelihood
if( any( a<low )|| any(a>high) )
{
like = -Inf;
} else
{
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
like = -Inf;
} else
{
#OJO
like = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
# Gradient
GG = cbind( filter(cbind(0, 1),cbind( 1, -a[ 3] ), matrix( 1, t, 1 )),
filter(cbind( 0, 1 ),cbind( 1, -a[3] ), y * (df-2) / df ),
filter(cbind( 0, 1 ),cbind( 1, -a[ 3 ]), h ));
colnames(GG)<-NULL;
g1 = ( ( df+1 ) * ( y / ( y + df * h ) ) - 1 ) / h;
#OJO
g1 = h / d + ( 2 + df ) * y / ( z + d * df) - ( 2 + df ) * h * z /( z + d * df ) / d;
G = GG * ( g1 %*% matrix( 1, 1, 3) );
gra = apply( G, 2, sum );
# Hessian
g2 = 1 / d + 2 * h^2 / d^2 - (2 + df) * y / (z + d * df)^2 * ( -2 * y - 2 * df * h) -
(2 + df) * z / ( z + d * df) / d + 2 * (2 + df) * h * y / ( z + d * df ) / d +
(2 + df) * h * z / (z + d * df)^2 / d *( -2 * y - 2 * df * h ) -
2 * ( 2 + df ) * h^2 * z / ( z + d * df ) / d ^2;
GG2 = GG[ ,c( 1, 2, 3, 1, 2, 3, 1, 2, 3)] * GG [, c( 1, 1, 1, 2, 2, 2, 3, 3, 3 ) ];
HH = matrix( 0, t, 9);
HH[ , 3 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ , 1 ] );
HH[ , 7 ] = HH[ , 3 ];
HH[ , 6 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ , 2 ] );
HH[ , 8 ] = HH[ ,6 ];
HH[ , 9 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ , 3 ] );
H = GG2 * (g2 %*% matrix( 1, 1, 9)) + HH * (g1 %*% matrix( 1, 1, 9)) ;
hes = matrix( apply( H, 2, sum), 3, 3 );
e = eigen(hes);
if( all(e$values>0) )
{
hes = -diag( 1, 3 );
negdef = 0;
}else
{
if( any(e$values > 0) )
{
negdef = 0;
val = pmin( e$values, max( e$values[ e$values<0 ] ) / big);
hes = e$vectors %*% diag( val, length(val)) %*% t(e$vectors);
} else
{
negdef = 1;
}
}
# Steepest Ascent or Newton
if( any( c(a==low, a==high) ) )
{
da = -((gra %*% t(gra))/(gra %*% hes %*% t(gra))) %*% gra;
}else
{
da = -gra %*% solve(hes);
}
# Termination criterion
term = (da %*% gra)[1] ;
if( ( term < tol1 ) && negdef )
{
break;
}
# If you are on the boundary and want to get out, slide along
if( length(da[( a==low ) & (da<0) ]))
{
da[( a==low ) & (da<0) ] = matrix( 0, dim(da[ ( a==low ) & ( da<0 ) ]));
}
if( length( da[( a==high ) & (da>0) ] ) )
{
da[( a==high ) & (da>0) ] = matrix( 0, dim(da[ ( a==high ) & ( da>0 ) ]));
}
# If you are stuck in a corner, terminate too
if( all(da==0) )
{
break;
}
# Go no further than next boundary
hit = cbind(( low[ da != 0 ] - a[ da != 0 ] ) / da[ da != 0 ],
( high[ da != 0 ] - a[ da != 0 ] ) / da[ da != 0 ]);
if( length(hit[hit <= 0])) hit[ hit <= 0 ] = 0;
da = min(cbind(hit, 1)) * da;
# Step search
best = 0;
newa = a + gold^( best - 1 ) * da;
if( any( newa < low ) || any(newa > high) )
{
left = -Inf;
}else
{
h = filter(cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) +
filter(cbind( 0, newa[ 1 ] ), cbind( 1, -newa[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
left = -Inf;
} else
{
left = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
newa = a + gold^( best ) * da;
if( any( newa < low ) || any(newa > high) )
{
center = -Inf;
}else
{
h = filter(cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) +
filter(cbind( 0, newa[ 1 ] ), cbind( 1, -newa[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
center = -Inf;
} else
{
center = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
newa = a + gold^( best + 1 ) * da;
if( any( newa < low ) || any(newa > high) )
{
right = -Inf;
}else
{
h = filter(cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) +
filter(cbind( 0, newa[ 1 ] ), cbind( 1, -newa[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
right = -Inf;
} else
{
right = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
if( all( like > c( left, center, right ) )|| all( left > c( center, right ) ))
{
while( 1 )
{
best = best - 1;
center = left;
newa = a + gold^( best - 1 ) * da;
if( any( newa < low ) || any(newa > high) )
{
left = -Inf;
}else
{
h = filter(cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) +
filter(cbind( 0, newa[ 1 ] ), cbind( 1, -newa[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
left = -Inf;
} else
{
left = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
if( all(center >= c( like, left ) ) )
{
break;
}
}
}else
{
if( all( right > c( left, center ) ) )
{
best = best + 1;
center = right;
newa = a + gold^( best + 1 ) * da;
while( 1 )
{
if( any( newa < low ) || any(newa > high) )
{
right = -Inf;
}else
{
h = filter(cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) +
filter(cbind( 0, newa[ 1 ] ), cbind( 1, -newa[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
right = -Inf;
} else
{
right = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
if( center > right )
{
break;
}
}
}
}
}
q = a;
return( q );
}
#
# @param A : [matrix] an indefinite symmetric matrix with non-negative diagonal elements
#
# @return XXX : [matrix] positive semi-definite matrix with same diagonal elements as A that is closest
# to A according to the Frobenius norm
#
# @note MATLAB's code written initially by Ilya Sharapov (1997)
#
# @references
# A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170}.
#
# See Meucci's script for "FitMultivariateGarch.m"
#
# @author Xavier Valls \email{flamejat@@gmail.com}
minfro = function( A )
{
if(!require("Matrix")) stop("Matrix package required for this script");
if( any( diag( A ) < 0) )
{
stop("Diagonal Elements Must Be Non-Negative!");
}else if( any( any( A != t(A) ) ) )
{
stop("Matrix Must Be Symmetric!");
}else if( all( eigen(A)$values >= 0 ) )
{
XXX = A;
}else
{
# if things go wrong make rho bigger and wait longer
rho = 0.75;
tol = 3e-6; # tolerance
maxj = 10; # max number of iterations
n = nrow( A );
M = diag( diag(A), ncol(A) ); # initialize with diagonal
oldnorm = norm( M - A , "F" );
oldnormj = oldnorm;
normj[ 1 ] = oldnorm;
j = 1;
incmax = 1e32; # just to enter the loop
while((j<maxj) && (incmax>tol))
{
incmax = 0;
for( i in 1 : n )
{
a = rbind( A[ 1:i-1, i ], A[ i+1:n, i ] );
m = rbind( M[ 1:i-1, i ], M[ i+1:n, i ] );
aii = A[ i, i ];
b = a - rho * m;
# Newton's step
x = newton( M, i, b, m, aii, n, rho );
P = as( diag( 1, n ), "sparseMatrix" );
P[ i, 1:n ] = t(x);
# update
Mtest = P %*% M %*% t(P);
M = Mtest;
inc = oldnorm - norm( M - A, "F" );
oldnorm = norm( M - A, "F" );
# find maximal increment over iteration
if( inc > incmax )
{
incmax = inc;
}
}
normj[ j+1 ] = oldnorm;
incj[ j ] = oldnormj - oldnorm;
oldnormj = oldnorm;
j = j + 1;
}
XXX = M;
}
return( XXX );
}
newton = function(M, i, b, m, aii, n, rho)
{
## Newton setp
# Subroutine called interbally by minfro.m
maxit = 40;
eps = 1e-9; # small correction
# to handle singularity
l = 0.0;
MM = rbind( cbind( M[ 1:i-1, 1:i-1 ], M[ 1:i-1, i+1:n ] ), cbind( M[ i+1:n, 1:i-1 ], M[ i+1:n, i+1:n ]) ) + eps * eye( n - 1 );
j = 1;
# loop
while( j < maxit )
{
tmp = MM %*% MM + l * MM;
IM = tmp %*% solve(diag( 1, dim( tmp )));
#IM = inv(MM*MM+l*MM);
x = IM %*% ( MM %*% b - l * rho * m);
f = rho * rho * aii + 2 * rho * t(x) %*% m + t(x) %*% MM %*% x - aii;
if( abs(f) < 1e-7 )
{
break;
}
dfdl = -2 * t( rho * m + MM %*% x ) %*% IM %*% ( rho * m + MM %*% x );
# Newton's step
l = l - f / dfdl;
j = j + 1;
}
if( abs(f) < 1e-7 )
{
# converged
xx = rbind( x[ 1:i-1 ], rho, x[ i:n-1 ]);
}else
{
# didn't converge
xx = matrix( 0, n, 1 );
xx[ i ] = 1;
}
return( xx );
}
R-Finance/Meucci documentation built on May 8, 2019, 3:52 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' Estimation of multivariate GARCH models
#'
#' @param returns : [matrix] (T x N) returns so rows must correspond to time and columns to assets
#' @param demean : [scalar] specifies whether returns should be demeaned (if demean = 1) or not to estimate the model; default value is 1.
#' @param eps : [scalar] used in enforcing a_ii + b_ii <= 1 - eps; the default value is zero
#' @param df : [scalar] degree of freedom for the t-distribution; the default value is 500 to make it, basically, normal
#'
#' @return mu : [vector]
#' @return ATMF : [matrix] coefficient matrix A-tilde (in the notation of the paper)
#' @return BTMF : [matrix] coefficient matrix B-tilde (in the notation of the paper)
#' @return CTMF : [matrix] coefficient matrix C-tilde (in the notation of the paper)
#' @return Hhat : [matrix] forecasted conditional covariance matrix
#'
#' @note Code for MATLAB initially written by Olivier Ledoit and Michael Wolf
#'
#' @references
#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170},
#' "E 136 - Equity market: multivariate GARCH process".
#'
#' See Meucci's script for "FitMultivariateGarch.m"
#'
#' @author Xavier Valls \email{flamejat@@gmail.com}
#' @export
FitMultivariateGarch = function( returns, demean = 1, eps = 0, df = 500 )
{
if (eps < 0)
{
error("eps must be a (small) positive number")
}
# Initialization
T = nrow( returns );
N = ncol( returns );
if ( 1 == demean )
{
mu = matrix(apply( returns, 2, mean ));
returns = returns-mu[ array(1,T ), ];
}
S = t(returns) %*% returns / ( T - 1 );
x = t(returns);
A = matrix( 0, N, N );
B = matrix( 0, N, N );
C = matrix( 0, N, N );
# Rescale Data
scale = sqrt( matrix( apply( t(x)^2, 2, mean ) ) );
x = x / scale[ , array( 1, T ) ];
# Estimation of On-Diagonal Elements
h = matrix( 0, N, T );
for( i in 1 : N )
{
# Likelihood Maximization
q = garch1f4( t( x[ i, ] ), eps, df)$q;
A[ i, i ] = q[ 2 ];
B[ i, i ] = q[ 3 ];
C[ i, i ] = q[ 1 ];
h[ i, ] = filter( cbind( 0, q[ 2] ), cbind( 1, -q[ 3 ] ), x[ i, ]^2 * (df-2) / df, mean( x[ i, ] ^ 2) * ( df-2 ) / df ) +
filter( cbind( 0, q[ 1 ] ), cbind( 1, -q[ 3 ] ), matrix( 1, 1, T) );
}
# First-step Estimation of Off-Diagonal Elements
for( i in 1 : (N-1) )
{
for( j in (i +1) : N )
{
# Likelihood Maximization
theta = garch2f8( x[ i, ] * x[ j, ], C[ i, i ], A[ i, i ], B[ i, i ], x[ i, ]^2, h[ i, ], C[ j, j ], A[ j, j ], B[ j, j ], x[ j, ]^2, h[ j, ], df );
A[ i, j ] = theta[ 2 ];
B[ i, j ] = theta[ 3 ];
C[ i, j ] = theta[ 1 ];
A[ j, i ] = A[ i, j ];
B[ j, i ] = B[ i, j ];
C[ j, i ] = C[ i, j ];
}
}
theta = garch2f8( x[ N, ] * x[ N, ], C[ N, N ], A[ N, N ], B[ N, N ], x[ N, ]^2, h[ N, ], C[ N, N ], A[ N, N ], B[ N, N ], x[ N, ]^2, h[ N, ], df );
A[ N, N ] = theta[ 2 ];
B[ N, N ] = theta[ 3 ];
C[ N, N ] = theta[ 1 ];
A[ N, N ] = A[ N, N ];
B[ N, N ] = B[ N, N ];
C[ N, N ] = C[ N, N ];
# Transformation of Coefficient Matrices
ATMF = minfro( A );
BTMF = minfro( B );
CTMF = minfro( C /( 1 - B )) * ( 1 - BTMF );
# Rescale
#C = C .* (scale * scale');
CTMF = CTMF * ( scale %*% t(scale) );
# Forecast of Conditional Covariance Matrix
Hhat = matrix( 0, N, N );
for( i in 1 : N )
{
for( j in 1 : N )
{
hSeries = filter( cbind( 0, ATMF[ i, j ] ), cbind( 1, -BTMF[ i, j ] ), returns[ , i ] * returns[ , j ], S[ i, j ] ) +
filter( cbind( 0, CTMF[ i, j ] ), cbind( 1, -BTMF[ i, j ] ), matrix( 1, T,1));
Hhat[ i, j ] = hSeries[ T ];
}
}
mu = matrix( mu );
return( list( mu = mu, ATMF = ATMF, BTMF = BTMF, CTMF = CTMF, Hhat = Hhat ) );
}
#' Fit a GARCH(1,1) model with student-t errors
#'
#' @param x : [vector] (T x 1) data generated by a GARCH(1,1) process
#' @param eps : [scalar] used in enforcing a_ii + b_ii <= 1 - eps; the default value is zero
#' @param df : [scalar] degree of freedom for the t-distribution; the default value is 500 to make it, basically, normal
#'
#' @return q : [vector] (4 x 1) parameters of the GARCH(1,1) process
#' @return qerr : [vector] (4 x 1) standard error of parameter estimates
#' @return hf : [scalar] current conditional heteroskedasticity estimate
#' @return hferr : [scalar] standard error on hf
#'
#' @note
#' MATLAB's script initially written by Olivier Ledoit, 4/28/1997
#'
#' Uses a conditional t-distribution with fixed degrees of freedom
#'
#' Difference with garch1f: errors come from the score alone
#'
#' @references
#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170}.
#'
#' See Meucci's script for "FitMultivariateGarch.m"
#'
#' @author Xavier Valls \email{flamejat@@gmail.com}
#' @export
garch1f4 = function( x, eps, df )
{
# Parameters
gold = ( 1 + sqrt(5) ) / 2; # step size increment
tol1 = 1e-7; # for termination criterion
tol2 = 1e-7; # for closeness to boundary
big = 2; # for making the hessian negative definite
maxiter = 50; # maximum number of iterations
n = 30; # number of points on the grid
# Rescale
y = ( array(x) - mean( array(x)))^2;
t = length(y);
scale = sqrt(mean( y^2 ));
y = y / scale;
s = mean(y);
# Grid search
#[ag, bg] = meshgrid( seq( 0, 1 - eps, length = n ), seq( 0, 1-eps, length = n ) );
ag = outer(seq( 0, 1 - eps, length = n )*0, seq( 0, 1-eps, length = n ), FUN = "+" );
bg = outer(seq( 0, 1 - eps, length = n ) , seq( 0, 1-eps, length = n ) * 0, FUN = "+" )
cg = pmax( s * (1-ag-bg), 0 );
likeg = matrix( -Inf, n, n );
for( i in 1 : n )
{
for( j in 1 : (n-i+1) )
{
h = filter( cbind( 0, ag[ i, j ] ), cbind( 1, -bg[ i, j] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, cg[ i, j ]),cbind( 1, - bg[ i, j ] ),matrix(1, t, 1));
likeg[ i, j ] = -sum( log(h) + (df+1) * log( 1 + y /h /df) );
}
}
maxlikeg = max(likeg, na.rm = TRUE);
maxima = which(likeg == maxlikeg); ##ok<MXFND>
maximum = max( maxima );
a = cbind(cg[ maximum ], ag[ maximum], bg[ maximum ]);
best = 0;
da = 0;
#term = 1;
#negdef = 0;
iter = 0;
while( iter < maxiter )
{
iter = iter + 1;
# New parameter1
a = a + gold^best * da;
# Conditional variance
h = filter(cbind( 0, a[ 2 ] ),cbind( 1, -a[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) + filter(cbind( 0, a[ 1 ] ),cbind( 1, -a[ 3 ] ), matrix(1, t, 1 ) );
# Likelihood
if( any( a<0 )||( (a[ 2 ] + a[ 3 ]) > (1 - eps) ))
{
like = -Inf;
} else
{
like = -sum( log( h )+( df + 1 ) * log( 1 + y /h /df) );
}
GG = cbind( filter(cbind(0, 1),cbind( 1, -a[ 3] ), matrix( 1, t, 1 )),
filter(cbind( 0, 1 ),cbind( 1, -a[3] ), y * (df-2) / df ),
filter(cbind( 0, 1 ),cbind( 1, -a[ 3 ]), h ));
colnames(GG)<-NULL;
g1 = ( ( df+1 ) * ( y / ( y + df * h ) ) - 1 ) / h;
G = GG * (g1 %*% matrix( 1, 1, 3));
gra = apply( G, 2, sum );
# Hessian
GG2 = GG[ ,c( 1, 2, 3, 1, 2, 3, 1, 2, 3)] * GG [, c( 1, 1, 1, 2, 2, 2, 3, 3, 3 ) ];
g2 = -((df+1) * ( y /( y + df*h) ) - 1 ) / h^2 - ( df * ( df + 1 ) )*(y / ( y+df * h ) ^2 / h );
HH = matrix( 0, t, 9);
HH[ , 3 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ ,1 ] );
HH[ , 7 ] = HH[ , 3 ];
HH[ , 6 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ , 2 ] );
HH[ , 8 ] = HH[ ,6 ];
HH[ , 9 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ , 3 ] );
H = GG2 * (g2 %*% matrix( 1, 1, 9)) + HH * (g1 %*% matrix( 1, 1, 9)) ;
hes = matrix( apply( H, 2, sum), 3, 3 );
e = eigen(hes);
if( any(e$values > 0) )
{
negdef = 0;
d = min( e$values, max( e$values[ e$values<0 ] ) / big);
hes = e$vectors %*% diag(d, length(d)) %*% t(e$vectors);
} else
{
negdef = 1;
}
# Direction
da = -gra %*% solve( hes );
# Termination criterion
term = (da %*% gra)[1];
if ((term < tol1) && negdef)
{
break;
}
best = 0;
newa = a + gold^(best-1) * da;
if( any(newa < 0 ) || ( newa[ 2 ] + newa[ 3 ] > (1-eps) ) )
{
left = -Inf;
} else
{
h = filter( cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, newa[ 1 ]),cbind( 1, - newa[ 3 ] ),matrix(1, t, 1));
left = -sum( log(h) + (df+1) * log( 1 + y /h /df) )
}
newa = a + gold^best * da;
if( any(newa < 0 ) || ( newa[ 2 ] + newa[ 3 ] > (1-eps) ) )
{
center = -Inf;
} else
{
h = filter( cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, newa[ 1 ]),cbind( 1, - newa[ 3 ] ),matrix(1, t, 1));
center = -sum( log(h) + (df+1) * log( 1 + y /h /df) )
}
newa = a + gold^( best+1 ) * da;
if( any(newa < 0 ) || ( newa[ 2 ] + newa[ 3 ] > (1-eps) ) )
{
right=-Inf;
} else
{
h = filter( cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, newa[ 1 ]),cbind( 1, - newa[ 3 ] ),matrix(1, t, 1));
right = -sum( log(h) + (df+1) * log( 1 + y /h /df) )
}
if( all(like > c( left, center, right)) || all( left > c( center, right )) )
{
while( 1 )
{
best = best-1;
center = left;
newa = a + gold^( best-1 ) * da;
if( any(newa < 0 ) || ( newa[ 2 ] + newa[ 3 ] > (1-eps) ) )
{
left = -Inf;
} else
{
h = filter( cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, newa[ 1 ]),cbind( 1, - newa[ 3 ] ),matrix(1, t, 1));
left = -sum( log(h) + (df+1) * log( 1 + y /h /df) )
}
if( all( center >= c( like, left ) ) )
{
break;
}
}
}else
{
if( all( right > c( left, center ) ) )
{
while( 1 )
{
best = best+1;
center = right;
newa = a + gold^(best+1) *da;
if( any(newa < 0 ) || ( newa[ 2 ] + newa[ 3 ] > (1-eps) ) )
{
right = -Inf;
} else
{
h = filter( cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * ( df - 2 )/ df ) + filter(cbind( 0, newa[ 1 ]),cbind( 1, - newa[ 3 ] ),matrix(1, t, 1));
right = -sum( log(h) + (df+1) * log( 1 + y /h /df) )
}
if( center > right )
{
break;
}
}
}
}
# If stuck at boundary then stop
if( ( center == like ) && ( any(a<tol2) || ( ( a[ 2 ] + a[ 3 ] ) > (1-tol2) ) ) )
{
break;
}
# end of optimization loop
}
if( length(a[ a < tol2 ]) )
{
a[ a < tol2 ] = matrix( 0, 1, length(a[ a<tol2 ]));
}
if( ( a[ 2 ] + a[ 3 ] )> ( 1-tol2 ) )
{
if( a[ 2 ]< ( 1 - tol2 ) )
{
a[ 2 ] = a[ 2 ] + (1 - a[ 2 ]-a[ 3 ] );
}else
{
a[ 3 ] = a[ 3 ] + ( 1 - a[ 2 ]- a[ 3 ] );
}
}
# Estimation error and volatility forecast
#aerr=solve(t(G)%*%G);
tmp = ( t(G) %*% G );
aerr = tmp %*% solve(diag( 1, dim(tmp)));
hf = a[ 1 ] + a[ 2 ] * y[ t ] * (df-2) / df + a[ 3 ] * h[ t ];
gf = c( 1, y[ t ], h[ t ] ) + a[ 3 ] * GG[ t, ];
hferr = gf %*% aerr %*% gf;
aerr = t( diag(aerr) );
# Revert to original scale
a[ 1 ] = a[ 1 ] * scale;
aerr[ 1 ] = aerr[ 1 ] * scale^2;
hf = hf* scale;
hferr = hferr * scale^2;
aerr = sqrt(aerr);
hferr = sqrt(hferr);
q = a;
qerr = aerr;
return( list( q = q, qerr = qerr, hf = hf, hferr = hferr ) );
}
#' Off-diagonal parameter estimation in bivariate GARCH(1,1) when diagonal parameters are given.
#'
#' @param y : [vector] (T x 1) data generated by a GARCH(1,1) process
#' @param c1 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix C
#' @param a1 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix A
#' @param b1 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix B
#' @param y1 : [vector] (T x 1) data generated by a GARCH(1,1) process
#' @param h1 : [vector] (T x 1) data generated by a GARCH(1,1) process
#' @param c2 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix C
#' @param a2 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix A
#' @param b2 : [scalar] diagonal parameter of the GARCH(1,1) process taken from matrix B
#' @param y2 : [vector] (T x 1) data generated by a GARCH(1,1) process
#' @param h2 : [vector] (T x 1) generated by a GARCH(1,1) process
#' @param df : [scalar] degree of freedom for the t-distribution; the default value is 500 to make it, basically, normal
#'
#' @return q : [vector] (4 x 1) parameters of the GARCH(1,1) process
#' @return qerr : [vector] (4 x 1) standard error of parameter estimates
#' @return hf : [scalar] current conditional heteroskedasticity estimate
#' @return hferr : [scalar] standard error on hf
#'
#' @note
#' MATLAB's code initially written by Olivier Ledoit, 4/28/1997
#'
#' Uses a conditional t-distribution with fixed degrees of freedom
#'
#' Steepest Ascent on boundary, Hessian off boundary, no grid search
#'
#' @references
#' A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170}.
#'
#' See Meucci's script for "FitMultivariateGarch.m"
#'
#' @author Xavier Valls \email{flamejat@@gmail.com}
#' @export
garch2f8 = function( y, c1, a1, b1, y1, h1, c2, a2, b2, y2, h2, df )
{
# Parameters
gold = ( 1 + sqrt( 5 ) ) / 2; # step size increment
tol1 = 1e-7; # for termination criterion
tol2 = 1e-7; # for closeness to boundary
big = 2; # for making the hessian negative definite
maxiter = 50; # maximum number of iterations
#n=30; # number of points on the grid
# Prepare
t = length( y );
y1 = array( y1 );
y2 = array( y2 );
y = array( y );
s = mean( y );
#s1 = mean(y1);
#s2 = mean(y2);
h1 = array( h1 );
h2 = array( h2 );
# Bounds
low = c(-sqrt(c1*c2), 0, 0 ) + tol2;
high = c( sqrt(c1*c2), sqrt(a1*a2), sqrt(b1*b2) ) - tol2;
# Starting Point
a0 = 0.9 * sqrt( a1 * a2 );
b0 = 0.9 * sqrt( b1 * b2 );
c0 = mean( y ) * ( 1 - a0 - b0 ) %*% ( df - 2 )/df;
c0 = sign( c0 ) * min( abs( c0 ), 0.9 * sqrt( c1*c2 ));
# Initialize optimization
a = cbind( c0, a0, b0 );
best = 0;
da = 0;
#term=1;
#negdef=0;
iter = 0;
while( iter < maxiter )
{
iter = iter + 1;
# New parameter1
a = a + gold^best * da;
# Conditional variance
h = filter(cbind( 0, a[ 2 ] ),cbind( 1, -a[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) + filter(cbind( 0, a[ 1 ] ),cbind( 1, -a[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
# Likelihood
if( any( a<low )|| any(a>high) )
{
like = -Inf;
} else
{
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
like = -Inf;
} else
{
#OJO
like = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
# Gradient
GG = cbind( filter(cbind(0, 1),cbind( 1, -a[ 3] ), matrix( 1, t, 1 )),
filter(cbind( 0, 1 ),cbind( 1, -a[3] ), y * (df-2) / df ),
filter(cbind( 0, 1 ),cbind( 1, -a[ 3 ]), h ));
colnames(GG)<-NULL;
g1 = ( ( df+1 ) * ( y / ( y + df * h ) ) - 1 ) / h;
#OJO
g1 = h / d + ( 2 + df ) * y / ( z + d * df) - ( 2 + df ) * h * z /( z + d * df ) / d;
G = GG * ( g1 %*% matrix( 1, 1, 3) );
gra = apply( G, 2, sum );
# Hessian
g2 = 1 / d + 2 * h^2 / d^2 - (2 + df) * y / (z + d * df)^2 * ( -2 * y - 2 * df * h) -
(2 + df) * z / ( z + d * df) / d + 2 * (2 + df) * h * y / ( z + d * df ) / d +
(2 + df) * h * z / (z + d * df)^2 / d *( -2 * y - 2 * df * h ) -
2 * ( 2 + df ) * h^2 * z / ( z + d * df ) / d ^2;
GG2 = GG[ ,c( 1, 2, 3, 1, 2, 3, 1, 2, 3)] * GG [, c( 1, 1, 1, 2, 2, 2, 3, 3, 3 ) ];
HH = matrix( 0, t, 9);
HH[ , 3 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ , 1 ] );
HH[ , 7 ] = HH[ , 3 ];
HH[ , 6 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ , 2 ] );
HH[ , 8 ] = HH[ ,6 ];
HH[ , 9 ] = filter(cbind(0, 1), cbind(1, -a[ 3 ] ),GG[ , 3 ] );
H = GG2 * (g2 %*% matrix( 1, 1, 9)) + HH * (g1 %*% matrix( 1, 1, 9)) ;
hes = matrix( apply( H, 2, sum), 3, 3 );
e = eigen(hes);
if( all(e$values>0) )
{
hes = -diag( 1, 3 );
negdef = 0;
}else
{
if( any(e$values > 0) )
{
negdef = 0;
val = pmin( e$values, max( e$values[ e$values<0 ] ) / big);
hes = e$vectors %*% diag( val, length(val)) %*% t(e$vectors);
} else
{
negdef = 1;
}
}
# Steepest Ascent or Newton
if( any( c(a==low, a==high) ) )
{
da = -((gra %*% t(gra))/(gra %*% hes %*% t(gra))) %*% gra;
}else
{
da = -gra %*% solve(hes);
}
# Termination criterion
term = (da %*% gra)[1] ;
if( ( term < tol1 ) && negdef )
{
break;
}
# If you are on the boundary and want to get out, slide along
if( length(da[( a==low ) & (da<0) ]))
{
da[( a==low ) & (da<0) ] = matrix( 0, dim(da[ ( a==low ) & ( da<0 ) ]));
}
if( length( da[( a==high ) & (da>0) ] ) )
{
da[( a==high ) & (da>0) ] = matrix( 0, dim(da[ ( a==high ) & ( da>0 ) ]));
}
# If you are stuck in a corner, terminate too
if( all(da==0) )
{
break;
}
# Go no further than next boundary
hit = cbind(( low[ da != 0 ] - a[ da != 0 ] ) / da[ da != 0 ],
( high[ da != 0 ] - a[ da != 0 ] ) / da[ da != 0 ]);
if( length(hit[hit <= 0])) hit[ hit <= 0 ] = 0;
da = min(cbind(hit, 1)) * da;
# Step search
best = 0;
newa = a + gold^( best - 1 ) * da;
if( any( newa < low ) || any(newa > high) )
{
left = -Inf;
}else
{
h = filter(cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) +
filter(cbind( 0, newa[ 1 ] ), cbind( 1, -newa[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
left = -Inf;
} else
{
left = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
newa = a + gold^( best ) * da;
if( any( newa < low ) || any(newa > high) )
{
center = -Inf;
}else
{
h = filter(cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) +
filter(cbind( 0, newa[ 1 ] ), cbind( 1, -newa[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
center = -Inf;
} else
{
center = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
newa = a + gold^( best + 1 ) * da;
if( any( newa < low ) || any(newa > high) )
{
right = -Inf;
}else
{
h = filter(cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) +
filter(cbind( 0, newa[ 1 ] ), cbind( 1, -newa[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
right = -Inf;
} else
{
right = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
if( all( like > c( left, center, right ) )|| all( left > c( center, right ) ))
{
while( 1 )
{
best = best - 1;
center = left;
newa = a + gold^( best - 1 ) * da;
if( any( newa < low ) || any(newa > high) )
{
left = -Inf;
}else
{
h = filter(cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) +
filter(cbind( 0, newa[ 1 ] ), cbind( 1, -newa[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
left = -Inf;
} else
{
left = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
if( all(center >= c( like, left ) ) )
{
break;
}
}
}else
{
if( all( right > c( left, center ) ) )
{
best = best + 1;
center = right;
newa = a + gold^( best + 1 ) * da;
while( 1 )
{
if( any( newa < low ) || any(newa > high) )
{
right = -Inf;
}else
{
h = filter(cbind( 0, newa[ 2 ] ), cbind( 1, -newa[ 3 ] ), y * ( df - 2 ) / df, s * (df-2) / df ) +
filter(cbind( 0, newa[ 1 ] ), cbind( 1, -newa[ 3 ] ), matrix(1, t, 1 ) );
d = h1 * h2 - h^2;
z = h2 * y1 + h1 * y2 - 2 * h * y;
if( any( d <= 0 ) || any( 1+ z / d / df <= 0 ) )
{
right = -Inf;
} else
{
right = -sum( log( d ) + ( 2 + df ) * log( 1 + z / d / df) /2 );
}
}
if( center > right )
{
break;
}
}
}
}
}
q = a;
return( q );
}
#
# @param A : [matrix] an indefinite symmetric matrix with non-negative diagonal elements
#
# @return XXX : [matrix] positive semi-definite matrix with same diagonal elements as A that is closest
# to A according to the Frobenius norm
#
# @note MATLAB's code written initially by Ilya Sharapov (1997)
#
# @references
# A. Meucci - "Exercises in Advanced Risk and Portfolio Management" \url{http://symmys.com/node/170}.
#
# See Meucci's script for "FitMultivariateGarch.m"
#
# @author Xavier Valls \email{flamejat@@gmail.com}
minfro = function( A )
{
if(!require("Matrix")) stop("Matrix package required for this script");
if( any( diag( A ) < 0) )
{
stop("Diagonal Elements Must Be Non-Negative!");
}else if( any( any( A != t(A) ) ) )
{
stop("Matrix Must Be Symmetric!");
}else if( all( eigen(A)$values >= 0 ) )
{
XXX = A;
}else
{
# if things go wrong make rho bigger and wait longer
rho = 0.75;
tol = 3e-6; # tolerance
maxj = 10; # max number of iterations
n = nrow( A );
M = diag( diag(A), ncol(A) ); # initialize with diagonal
oldnorm = norm( M - A , "F" );
oldnormj = oldnorm;
normj[ 1 ] = oldnorm;
j = 1;
incmax = 1e32; # just to enter the loop
while((j<maxj) && (incmax>tol))
{
incmax = 0;
for( i in 1 : n )
{
a = rbind( A[ 1:i-1, i ], A[ i+1:n, i ] );
m = rbind( M[ 1:i-1, i ], M[ i+1:n, i ] );
aii = A[ i, i ];
b = a - rho * m;
# Newton's step
x = newton( M, i, b, m, aii, n, rho );
P = as( diag( 1, n ), "sparseMatrix" );
P[ i, 1:n ] = t(x);
# update
Mtest = P %*% M %*% t(P);
M = Mtest;
inc = oldnorm - norm( M - A, "F" );
oldnorm = norm( M - A, "F" );
# find maximal increment over iteration
if( inc > incmax )
{
incmax = inc;
}
}
normj[ j+1 ] = oldnorm;
incj[ j ] = oldnormj - oldnorm;
oldnormj = oldnorm;
j = j + 1;
}
XXX = M;
}
return( XXX );
}
newton = function(M, i, b, m, aii, n, rho)
{
## Newton setp
# Subroutine called interbally by minfro.m
maxit = 40;
eps = 1e-9; # small correction
# to handle singularity
l = 0.0;
MM = rbind( cbind( M[ 1:i-1, 1:i-1 ], M[ 1:i-1, i+1:n ] ), cbind( M[ i+1:n, 1:i-1 ], M[ i+1:n, i+1:n ]) ) + eps * eye( n - 1 );
j = 1;
# loop
while( j < maxit )
{
tmp = MM %*% MM + l * MM;
IM = tmp %*% solve(diag( 1, dim( tmp )));
#IM = inv(MM*MM+l*MM);
x = IM %*% ( MM %*% b - l * rho * m);
f = rho * rho * aii + 2 * rho * t(x) %*% m + t(x) %*% MM %*% x - aii;
if( abs(f) < 1e-7 )
{
break;
}
dfdl = -2 * t( rho * m + MM %*% x ) %*% IM %*% ( rho * m + MM %*% x );
# Newton's step
l = l - f / dfdl;
j = j + 1;
}
if( abs(f) < 1e-7 )
{
# converged
xx = rbind( x[ 1:i-1 ], rho, x[ i:n-1 ]);
}else
{
# didn't converge
xx = matrix( 0, n, 1 );
xx[ i ] = 1;
}
return( xx );
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.