一、
garchFit
函數的參數---------------------------------------------
algorithm
a string parameter that determines the algorithm used for maximum likelihood estimation.
設定最大似然估計所用的算法
cond.dist
a character string naming the desired conditional distribution. Valid values are "dnorm", "dged", "dstd", "dsnorm", "dsged", "dsstd" and "QMLE". The default value is the normal distribution. See Details for more information.
條件擾動(新息)的分布
control
control parameters, the same as used for the functions from nlminb, and 'bfgs' and 'Nelder-Mead' from optim.
data
an optional timeSeries or data frame object containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which armaFit is called. If data is an univariate series, then the series is converted into a numeric vector and the name of the response in the formula will be neglected(被忽略).
包含變量的
timeSeries或者數據框類型，如果打他中沒有這些變量，則將去環境(
formula
)中找,
armaFit的環境通常也會被調用。如果是單變量序列，則序列被轉為數值向量，
formula那些就被忽略。
delta
a numeric value, the exponent delta of the variance recursion. By default, this value will be fixed, otherwise the exponent will be estimated together with the other model parameters if include.delta=FALSE.
變量定義的delta指數？
description
a character string which allows for a brief description.
簡短描述
formula
formula object describing the mean and variance equation of the ARMA-GARCH/APARCH model. A pure GARCH(1,1) model is selected when e.g. formula=~garch(1,1). To specify for example an ARMA(2,1)-APARCH(1,1) use formula = ~arma(2,1)+apaarch(1,1).
描述
ARMA-GARCH/APARCH模型的均值方程和波動率方程的對象
~garch(1,1) ——表示
GARCH(1,1)
模型
~arma(2,1)+apaarch(1,1) ——表示
ARMA(2,1)-APARCH(1,1)模型
gamma
APARCH leverage parameter entering into the formula for calculating the expectation value.
杠桿參數
hessian
a string denoting how the Hessian matrix should be evaluated, either hessian ="rcd", or "ropt", the default, "rcd" is a central difference approximation implemented in R and "ropt" use the internal R function optimhess.
海賽矩陣被計算的方式，默認為
rcd(中心差分近似？)
include.delta
a logical flag which determines if the parameter for the recursion equation delta will be estimated or not. If include.delta=FALSE then the shape parameter will be kept fixed during the process of parameter optimization.
邏輯標志，表示迭代法處delta參數是否被估計
include.mean
this flag determines if the parameter for the mean will be estimated or not. If include.mean=TRUE this will be the case, otherwise the parameter will be kept fixed durcing(應該是during) the process of parameter optimization.
include.mean=TRUE表示均值是否被估計，否則均值參數在參數最優化過程中將是固定的
include.shape
a logical flag which determines if the parameter for the shape of the conditional distribution will be estimated or not. If include.shape=FALSE then the shape parameter will be kept fixed during the process of parameter optimization.
一個參數，決定條件分布的自由度shape是否被估計
include.shape=FALSE則在參數最優化過程中，
t分布的自由度shape是固定的
include.skew
a logical flag which determines if the parameter for the skewness of the conditional distribution will be estimated or not. If include.skew=FALSE then the skewness parameter will be kept fixed during the process of parameter optimization.
skew=FALSE則shew(偏度)在參數最優化過程中是固定的
init.rec
a character string indicating the method how to initialize the mean and varaince recursion relation.
一個字符串，指明如何初始化均值和方差迭代關系。
取值呢？
"mci",
"uev"
leverage
a logical flag for APARCH models. Should the model be leveraged? By default leverage=TRUE.
是否帶杠桿
shape
a numeric value, the shape parameter of the conditional distribution.
條件分布的自由度值
skew
a numeric value, the skewness parameter of the conditional distribution.
條件分布的偏度
title
a character string which allows for a project title.
標題
trace
a logical flag. Should the optimization process of fitting the model parameters be printed? By default trace=TRUE.
參數最優化過程是否被打印出來。
...
additional arguments to be passed.
二、查看數據----------------------------------------------
雙擊da(相當於View(da))，可以看到數據是這樣的
經過對數轉化后的收益率，是我們建模用的原始數據。
intc=log(da$intc+1)intc [1] 0.0099998346 -0.1500127529 0.0670640792 [4] 0.0829486352 -0.1103484905 0.1251628488 [7] 0.4855078158 0.1112255825 0.2109235908..........
library(fGarch)da=read.table("m-intcsp7309.txt",header=T)intc=log(da$intc+1)m4=garchFit(~1+garch(1,1),data=intc)summary(m4)Title: GARCH Modelling Call: garchFit(formula = ~1 + garch(1, 1), data = intc) Mean and Variance Equation: data ~ 1 + garch(1, 1)<environment: 0x000000001d475658> [data = intc]Conditional Distribution: norm Coefficient(s): mu omega alpha1 beta1 0.01126568 0.00091902 0.08643831 0.85258554 Std. Errors: based on Hessian Error Analysis: Estimate Std. Error t value Pr(>|t|) mu 0.0112657 0.0053931 2.089 0.03672 * omega 0.0009190 0.0003888 2.364 0.01808 * alpha1 0.0864383 0.0265439 3.256 0.00113 ** beta1 0.8525855 0.0394322 21.622 < 2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Log Likelihood: 312.3307 normalized: 0.7034475 Description: Fri Jan 27 21:49:06 2017 by user: xuan Standardised Residuals Tests: Statistic p-Value Jarque-Bera Test R Chi^2 174.904 0 Shapiro-Wilk Test R W 0.9709615 1.030282e-07 Ljung-Box Test R Q(10) 8.016844 0.6271916 Ljung-Box Test R Q(15) 15.5006 0.4159946 Ljung-Box Test R Q(20) 16.41549 0.6905368 Ljung-Box Test R^2 Q(10) 0.8746345 0.9999072 Ljung-Box Test R^2 Q(15) 11.35935 0.7267295 Ljung-Box Test R^2 Q(20) 12.55994 0.8954573 LM Arch Test R TR^2 10.51401 0.5709617 Information Criterion Statistics: AIC BIC SIC HQIC -1.388877 -1.351978 -1.389037 -1.374326
三、主要的函數結構
----------------------------------------------
garchFit函數
match.xxx函數....
各種判斷可以忽略
.garchArgsParser
.garchFit函數
ans最終的返回結果
.garchFit函數
.garchInitSeries函數
.garchInitParameters函數
.garchOptimizeLLH函數
.garchRnlminb函數
nlminb函數
.garchLLH函數
.Fortran("garchllh")
四、調試
----------------------------------------------
(4.0)garchFit函數
function (formula = ~garch(1, 1), data = dem2gbp, init.rec = c("mci", "uev"), delta = 2, skew = 1, shape = 4, cond.dist = c("norm", "snorm", "ged", "sged", "std", "sstd", "snig", "QMLE"), include.mean = TRUE, include.delta = NULL, include.skew = NULL, include.shape = NULL, leverage = NULL, trace = TRUE, algorithm = c("nlminb", "lbfgsb", "nlminb+nm", "lbfgsb+nm"), hessian = c("ropt", "rcd"), control = list(), title = NULL, description = NULL, ...) { DEBUG = FALSE init.rec = match.arg(init.rec) cond.dist = match.arg(cond.dist) hessian = match.arg(hessian) algorithm = match.arg(algorithm) CALL = match.call() Name = capture.output(substitute(data)) if (is.character(data)) { eval(parse(text = paste("data(", data, ")"))) data = eval(parse(text = data)) } data <- as.data.frame(data) if (isUnivariate(data)) { colnames(data) <- "data" } else { uniqueNames = unique(sort(colnames(data))) if (is.null(colnames(data))) { stop("Column names of data are missing.") } if (length(colnames(data)) != length(uniqueNames)) { stop("Column names of data are not unique.") } } if (length(formula) == 3 && isUnivariate(data)) formula[2] <- NULL if (length(formula) == 2) { if (isUnivariate(data)) { formula = as.formula(paste("data", paste(formula, collapse = " "))) } else { stop("Multivariate data inputs require lhs for the formula.") } } robust.cvar <- (cond.dist == "QMLE") args = .garchArgsParser(formula = formula, data = data, trace = FALSE) if (DEBUG) print(list(formula.mean = args$formula.mean, formula.var = args$formula.var, series = args$series, init.rec = init.rec, delta = delta, skew = skew, shape = shape, cond.dist = cond.dist, include.mean = include.mean, include.delta = include.delta, include.skew = include.skew, include.shape = include.shape, leverage = leverage, trace = trace, algorithm = algorithm, hessian = hessian, robust.cvar = robust.cvar, control = control, title = title, description = description)) ans = .garchFit(formula.mean = args$formula.mean, formula.var = args$formula.var, series = args$series, init.rec, delta, skew, shape, cond.dist, include.mean, include.delta, include.skew, include.shape, leverage, trace, algorithm, hessian, robust.cvar, control, title, description, ...) ans@call = CALL attr(formula, "data") <- paste("data = ", Name, sep = "") ans@formula = formula ans}
以下四句都是對傳入的參數是否math(合法、規范)做檢驗的
init.rec = match.arg(init.rec) cond.dist = match.arg(cond.dist) hessian = match.arg(hessian) algorithm = match.arg(algorithm) CALL = match.call()
我們可以按下
進入
函數查看
進入
函數查看
截取match.arg函數中的幾句作說明：
formal.args <- formals(sys.function(sys.parent())).......................................if (identical(arg, choices)) return(arg[1L])
如果傳入和默認的一樣，說明沒有傳入，就返回第一個作為默認值
可以點擊
跳出函數
跳出函數
(4.0.1)
.garchArgsParser函數
args = .garchArgsParser(formula = formula, data = data, trace = FALSE)
function (formula, data, trace = FALSE) { allVars = unique(sort(all.vars(formula))) allVarsTest = mean(allVars %in% colnames(data)) if (allVarsTest != 1) { print(allVars) print(colnames(data)) stop("Formula and data units do not match.") } formula.lhs = as.character(formula)[2] mf = match.call(expand.dots = FALSE) if (trace) { cat("\nMatched Function Call:\n ") print(mf) } m = match(c("formula", "data"), names(mf), 0) mf = mf[c(1, m)] mf[[1]] = as.name(".garchModelSeries") mf$fake = FALSE mf$lhs = TRUE if (trace) { cat("\nModelSeries Call:\n ") print(mf) } x = eval(mf, parent.frame()) if (trace) print(x) x = as.vector(x[, 1]) names(x) = rownames(data) if (trace) print(x) allLabels = attr(terms(formula), "term.labels") if (trace) { cat("\nAll Term Labels:\n ") print(allLabels) } if (length(allLabels) == 2) { formula.mean = as.formula(paste("~", allLabels[1])) formula.var = as.formula(paste("~", allLabels[2])) } else if (length(allLabels) == 1) { formula.mean = as.formula("~ arma(0, 0)") formula.var = as.formula(paste("~", allLabels[1])) } if (trace) { cat("\nMean Formula:\n ") print(formula.mean) cat("\nVariance Formula:\n ") print(formula.var) } ans <- list(formula.mean = formula.mean, formula.var = formula.var, formula.lhs = formula.lhs, series = x) ans}
提下這個函數里面的一些信息:
data是數據
%in%運算符就是math的作用 ?'
%in%
'
x就是原始數據
(需要特別說明的是,x是從data里取出來后as.vector的，被轉化為向量后，會四舍五入顯示了，所以如果直接在environment pane中查看有點不橡原始數據，比如intc的第一個值是
0.0099998346
，而這里x的值顯示是0.01)
函數運行到最后一行的時候，
environment pane中的部分信息。
(方程的基本形式已經被解析出來了)
(4.0.2)
.garchFit函數
ans = .garchFit(formula.mean = args$formula.mean, formula.var = args$formula.var, series = args$series, init.rec, delta, skew, shape, cond.dist, include.mean, include.delta, include.skew, include.shape, leverage, trace, algorithm, hessian, robust.cvar, control, title, description, ...)
我們為什么分析了
.garchArgsParser
這個解析函數，就是因為
.garchFit
函數用到的許多參數都是從args中取出的。
另外提到一下，帶.開頭的函數和變量都是內部的，不會直接顯示出來的。
如果是變量名.xxx，這樣的形式，那么需要我們自己在控制台輸入
.xxx查看，不能在envir pane中查看。
如果要查看內部函數，則
單擊可以切換到對應的源代碼和變量
可以選擇環境，一般情況是在global environment
.garchFit函數的源碼
function (formula.mean = ~arma(0, 0), formula.var = ~garch(1, 1), series, init.rec = c("mci", "uev"), delta = 2, skew = 1, shape = 4, cond.dist = c("norm", "snorm", "ged", "sged", "std", "sstd", "QMLE"), include.mean = TRUE, include.delta = NULL, include.skew = NULL, include.shape = NULL, leverage = NULL, trace = TRUE, algorithm = c("sqp", "nlminb", "lbfgsb", "nlminb+nm", "lbfgsb+nm"), hessian = c("ropt", "rcd"), robust.cvar, control = list(), title = NULL, description = NULL, ...) { DEBUG <- FALSE if (DEBUG) print("Formula Specification ...") fcheck = rev(all.names(formula.mean))[1] if (fcheck == "ma") { stop("Use full formula: arma(0,q) for ma(q)") } else if (fcheck == "ar") { stop("Use full formula expression: arma(p,0) for ar(p)") } if (DEBUG) print("Recursion Initialization ...") if (init.rec[1] != "mci" & algorithm[1] != "sqp") { stop("Algorithm only supported for mci Recursion") } .StartFit <- Sys.time() if (DEBUG) print("Generate Control List ...") con <- .garchOptimizerControl(algorithm, cond.dist) con[(namc <- names(control))] <- control if (DEBUG) print("Initialize Time Series ...") data <- series scale <- if (con$xscale) sd(series) else 1 series <- series/scale .series <- .garchInitSeries(formula.mean = formula.mean, formula.var = formula.var, cond.dist = cond.dist[1], series = series, scale = scale, init.rec = init.rec[1], h.start = NULL, llh.start = NULL, trace = trace) .setfGarchEnv(.series = .series) if (DEBUG) print("Initialize Model Parameters ...") .params <- .garchInitParameters(formula.mean = formula.mean, formula.var = formula.var, delta = delta, skew = skew, shape = shape, cond.dist = cond.dist[1], include.mean = include.mean, include.delta = include.delta, include.skew = include.skew, include.shape = include.shape, leverage = leverage, algorithm = algorithm[1], control = con, trace = trace) .setfGarchEnv(.params = .params) if (DEBUG) print("Select Conditional Distribution ...") .setfGarchEnv(.garchDist = .garchSetCondDist(cond.dist[1])) if (DEBUG) print("Estimate Model Parameters ...") .setfGarchEnv(.llh = 1e+99) .llh <- .getfGarchEnv(".llh") fit = .garchOptimizeLLH(hessian, robust.cvar, trace) if (DEBUG) print("Add to fit ...") .series <- .getfGarchEnv(".series") .params <- .getfGarchEnv(".params") names(.series$h) <- NULL fit$series = .series fit$params = .params if (DEBUG) print("Retrieve Residuals and Fitted Values ...") residuals = .series$z fitted.values = .series$x - residuals h.t = .series$h if (.params$includes["delta"]) deltainv = 1/fit$par["delta"] else deltainv = 1/fit$params$delta sigma.t = (.series$h)^deltainv if (DEBUG) print("Standard Errors and t-Values ...") fit$cvar <- if (robust.cvar) (solve(fit$hessian) %*% (t(fit$gradient) %*% fit$gradient) %*% solve(fit$hessian)) else -solve(fit$hessian) fit$se.coef = sqrt(diag(fit$cvar)) fit$tval = fit$coef/fit$se.coef fit$matcoef = cbind(fit$coef, fit$se.coef, fit$tval, 2 * (1 - pnorm(abs(fit$tval)))) dimnames(fit$matcoef) = list(names(fit$tval), c(" Estimate", " Std. Error", " t value", "Pr(>|t|)")) if (DEBUG) print("Add Title and Description ...") if (is.null(title)) title = "GARCH Modelling" if (is.null(description)) description = description() Time = Sys.time() - .StartFit if (trace) { cat("\nTime to Estimate Parameters:\n ") print(Time) } new("fGARCH", call = as.call(match.call()), formula = as.formula(paste("~", formula.mean, "+", formula.var)), method = "Max Log-Likelihood Estimation", data = data, fit = fit, residuals = residuals, fitted = fitted.values, h.t = h.t, sigma.t = as.vector(sigma.t), title = as.character(title), description = as.character(description))}
在正式分析函數的源碼之前
先看一下模型變量m4的結構：
> str(m4)
Formal class 'fGARCH' [package "fGarch"] with 11 slots
..@ call : language garchFit(formula = ~1 + garch(1, 1), data = intc, trace = F)
..@ formula :Class 'formula' language data ~ 1 + garch(1, 1)
.. .. ..- attr(*, ".Environment")=<environment: 0x0000000003e60360>
.. .. ..- attr(*, "data")= chr "data = intc"
..@ method : chr "Max Log-Likelihood Estimation"
..@ data : Named num [1:444] 0.01 -0.15 0.0671 0.0829 -0.1103 ...
.. ..- attr(*, "names")= chr [1:444] "1" "2" "3" "4" ...
..@ fit :List of 17
.. ..$ par : Named num [1:4] 0.011266 0.000919 0.086438 0.852586 #參數
.. .. ..- attr(*, "names")= chr [1:4] "mu" "omega" "alpha1" "beta1"
.. ..$ objective : num 604
.. ..$ convergence: int 1
.. ..$ iterations : int 19
.. ..$ evaluations: Named int [1:2] 36 105
.. .. ..- attr(*, "names")= chr [1:2] "function" "gradient"
.. ..$ message : chr "singular convergence (7)"
.. ..$ value : Named num -312 (為什么顯示是負的，而輸出是正的？，因為nlminb是取最小值的)
.. .. ..- attr(*, "names")= chr "LogLikelihood"
.. ..$ coef : Named num [1:4] 0.011266 0.000919 0.086438 0.852586
.. .. ..- attr(*, "names")= chr [1:4] "mu" "omega" "alpha1" "beta1"
.. ..$ llh : Named num -312 #似然函數的最大值取負數
.. .. ..- attr(*, "names")= chr "LogLikelihood"
.. ..$ hessian : num [1:4, 1:4] -35115 108140 -253 919 108140 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:4] "mu" "omega" "alpha1" "beta1"
.. .. .. ..$ : chr [1:4] "mu" "omega" "alpha1" "beta1"
.. .. ..- attr(*, "time")=Class 'difftime' atomic [1:1] 0.021
.. .. .. .. ..- attr(*, "units")= chr "secs"
.. ..$ ics : Named num [1:4] -1.39 -1.35 -1.39 -1.37
.. .. ..- attr(*, "names")= chr [1:4] "AIC" "BIC" "SIC" "HQIC"
.. ..$ series :List of 10
.. .. ..$ model : chr [1:2] "arma" "garch"
.. .. ..$ order : Named num [1:4] 0 0 1 1
.. .. .. ..- attr(*, "names")= chr [1:4] "u" "v" "p" "q"
.. .. ..$ max.order: num 1
.. .. ..$ z : num [1:444] -0.00127 -0.16128 0.0558 0.07168 -0.12161 ... 殘差
.. .. ..$ h : num [1:444] 0.016 0.0146 0.0156 0.0145 0.0137 ...
.. .. ..$ x : Named num [1:444] 0.01 -0.15 0.0671 0.0829 -0.1103 ...
.. .. .. ..- attr(*, "names")= chr [1:444] "1" "2" "3" "4" ...
.. .. ..$ scale : num 0.127
.. .. ..$ init.rec : chr "mci"
.. .. ..$ h.start : num 2
.. .. ..$ llh.start: num 1
.. ..$ params :List of 14
.. .. ..$ params : Named num [1:8] 0.011266 0.000919 0.086438 0.1 0.852586 ...
.. .. .. ..- attr(*, "names")= chr [1:8] "mu" "omega" "alpha1" "gamma1" ...
.. .. ..$ U : Named num [1:8] -1.13 1.00e-06 1.00e-08 -1.00 1.00e-08 ...
.. .. .. ..- attr(*, "names")= chr [1:8] "mu" "omega" "alpha1" "gamma1" ...
.. .. ..$ V : Named num [1:8] 1.13 100 1 1 1 ...
.. .. .. ..- attr(*, "names")= chr [1:8] "mu" "omega" "alpha1" "gamma1" ...
.. .. ..$ includes : Named logi [1:8] TRUE TRUE TRUE FALSE TRUE FALSE ...
.. .. .. ..- attr(*, "names")= chr [1:8] "mu" "omega" "alpha1" "gamma1" ...
.. .. ..$ index : Named int [1:4] 1 2 3 5
.. .. .. ..- attr(*, "names")= chr [1:4] "mu" "omega" "alpha1" "beta1"
.. .. ..$ mu : Named num 0.113
.. .. .. ..- attr(*, "names")= chr "mu"
.. .. ..$ delta : num 2 #默認是2
.. .. ..$ skew : num 1
.. .. ..$ shape : num 4
.. .. ..$ cond.dist : chr "norm"
.. .. ..$ leverage : logi FALSE
.. .. ..$ persistence: Named num 0.9
.. .. .. ..- attr(*, "names")= chr "persistence"
.. .. ..$ control :List of 19
.. .. .. ..$ fscale : logi TRUE
.. .. .. ..$ xscale : logi TRUE
.. .. .. ..$ algorithm: chr "nlminb"
.. .. .. ..$ llh : chr "internal"
.. .. .. ..$ tol1 : num 1
.. .. .. ..$ tol2 : num 1
.. .. .. ..$ MIT : num 2000
.. .. .. ..$ MFV : num 5000
.. .. .. ..$ MET : num 5
.. .. .. ..$ MEC : num 2
.. .. .. ..$ MER : num 1
.. .. .. ..$ MES : num 4
.. .. .. ..$ XMAX : num 1000
.. .. .. ..$ TOLX : num 1e-10
.. .. .. ..$ TOLC : num 1e-06
.. .. .. ..$ TOLG : num 1e-06
.. .. .. ..$ TOLD : num 1e-06
.. .. .. ..$ TOLS : num 1e-04
.. .. .. ..$ RPF : num 0.01
.. .. ..$ llh : num 604
.. ..$ cvar : num [1:4, 1:4] 2.91e-05 1.06e-07 -1.51e-05 6.60e-06 1.06e-07 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:4] "mu" "omega" "alpha1" "beta1"
.. .. .. ..$ : chr [1:4] "mu" "omega" "alpha1" "beta1"
.. ..$ se.coef : Named num [1:4] 0.005393 0.000389 0.026544 0.039432
.. .. ..- attr(*, "names")= chr [1:4] "mu" "omega" "alpha1" "beta1"
.. ..$ tval : Named num [1:4] 2.09 2.36 3.26 21.62
.. .. ..- attr(*, "names")= chr [1:4] "mu" "omega" "alpha1" "beta1"
.. ..$ matcoef : num [1:4, 1:4] 0.011266 0.000919 0.086438 0.852586 0.005393 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:4] "mu" "omega" "alpha1" "beta1"
.. .. .. ..$ : chr [1:4] " Estimate" " Std. Error" " t value" "Pr(>|t|)"
..@ residuals : num [1:444] -0.00127 -0.16128 0.0558 0.07168 -0.12161 ...
..@ fitted : Named num [1:444] 0.0113 0.0113 0.0113 0.0113 0.0113 ...
#是
.series$x - residuals 其實這個x就是r(i)，
#那么fitted就是mu了(只不過這里的值是被四舍五入處理了不明顯)
代碼
fitted.values = .series$x - residuals
fitted = fitted.values
.. ..- attr(*, "names")= chr [1:444] "1" "2" "3" "4" ...
..@ h.t : num [1:444] 0.016 0.0146 0.0156 0.0145 0.0137 ... #條件方差
..@ sigma.t : num [1:444] 0.127 0.121 0.125 0.12 0.117 ... #是條件標准差，等於volatility(m4)
代碼
deltainv = 1/fit$par["delta"]=0.
sigma.t = (.series$h)^deltainv
..@ title : chr "GARCH Modelling"
..@ description: chr "Sun Jan 22 10:09:46 2017 by user: xuan"
?'@'
Extract or replace the contents of a slot in a object with a formal (S4) class structure.
從S4類里取出slot可以適應@運算符
hessian
a string denoting how the Hessian matrix should be evaluated, either hessian ="rcd", or "ropt", the default, "rcd" is a central difference approximation implemented in R and "ropt" use the internal R function optimhess.
黑塞矩陣（Hessian Matrix），又譯作海森矩陣、海瑟矩陣、海塞矩陣等，是一個多元函數的二階偏導數構成的方陣，描述了函數的局部曲率
Details
"QMLE" stands for Quasi-Maximum Likelihood Estimation, which assumes normal distribution and uses robust standard errors for inference.
Bollerslev and Wooldridge (1992) proved that if the mean and the volatility equations are correctly specified, the QML estimates are consistent and asymptotically normally distributed. However, the estimates are not efficient and “the efficiency loss can be marked under asymmetric ... distributions” (Bollerslev and Wooldridge (1992), p. 166).
The robust variance-covariance matrix of the estimates equals the (Eicker-White) sandwich estimator, i.e.
V = H^(-1) G' G H^(-1),
where V denotes the variance-covariance matrix, H stands for the Hessian and G represents the matrix of contributions to the gradient(傾斜度), the elements of which are defined as
G_{t,i} = derivative of l_{t} w.r.t. zeta_{i},
where l_{t} is the log likelihood of the t-th observation and zeta_{i} is the i-th estimated parameter. See sections 10.3 and 10.4 in Davidson and MacKinnon (2004) for a more detailed description of the robust variance-covariance matrix.
Value
garchFit
returns a S4 object of class "fGARCH" with the following slots:
@call
the call of the garch function.
@formula
a list with two formula entries, one for the mean and the other one for the variance equation.
@method
a string denoting the optimization method, by default the returneds string is "Max Log-Likelihood Estimation".
@data
a list with one entry named x, containing the data of the time series to be estimated, the same as given by the input argument series.
@fit
a list with the results from the parameter estimation. The entries of the list depend on the selected algorithm, see below.
@residuals
a numeric vector with the (raw, unstandardized) residual values.
@fitted
a numeric vector with the fitted values.
@h.t
a numeric vector with the conditional variances (h.t = sigma.t^delta).
@sigma.t
a numeric vector with the conditional standard deviation.
@title
a title string.
@description
a string with a brief description.
The entries of the @fit slot show the results from the optimization.
其中1
scale <- if (con$xscale) sd(series)series <- series/scale
是對數據進行標准化
（注意，標准化之后，數據和原始數據不同啦，在后面運算時，不要忘了~）
我們可以在調試界面的控制台直接輸入命令，來查看變量在當前環境和一定語句代運行完后的值
這個所謂的標准化，並沒有減去均值的，要注意！
其中2
.series <- .garchInitSeries(formula.mean = formula.mean, formula.var = formula.var, cond.dist = cond.dist[1], series = series, scale = scale, init.rec = init.rec[1], h.start = NULL, llh.start = NULL, trace = trace).setfGarchEnv(.series = .series)
注意這種.setfGarchEnv的做法，就是對環境變量中的變量賦值
還有.getfGarchEnv，可以提高自己的編程能力。
.garchInitSeries函數
function (formula.mean, formula.var, cond.dist, series, scale, init.rec, h.start, llh.start, trace) { mm = length(formula.mean) if (mm != 2) stop("Mean Formula misspecified") end = regexpr("\\(", as.character(formula.mean[mm])) - 1 model.mean = substr(as.character(formula.mean[mm]), 1, end) if (!any(c("ar", "ma", "arma") == model.mean)) stop("formula.mean must be one of: ar, ma, arma") mv = length(formula.var) if (mv != 2) stop("Variance Formula misspecified") end = regexpr("\\(", as.character(formula.var[mv])) - 1 model.var = substr(as.character(formula.var[mv]), 1, end) if (!any(c("garch", "aparch") == model.var)) stop("formula.var must be one of: garch, aparch") model.order = as.numeric(strsplit(strsplit(strsplit(as.character(formula.mean), "\\(")[[2]][2], "\\)")[[1]], ",")[[1]]) u = model.order[1] v = 0 if (length(model.order) == 2) v = model.order[2] maxuv = max(u, v) if (u < 0 | v < 0) stop("*** ARMA orders must be positive.") model.order = as.numeric(strsplit(strsplit(strsplit(as.character(formula.var), "\\(")[[2]][2], "\\)")[[1]], ",")[[1]]) p = model.order[1] q = 0 if (length(model.order) == 2) q = model.order[2] if (p + q == 0) stop("Misspecified GARCH Model: Both Orders are zero!") maxpq = max(p, q) if (p < 0 | q < 0) stop("*** GARCH orders must be positive.") max.order = max(maxuv, maxpq) if (is.null(h.start)) h.start = max.order + 1 if (is.null(llh.start)) llh.start = 1 if (init.rec != "mci" & model.var != "garch") { stop("Algorithm only supported for mci Recursion") } ans = list(model = c(model.mean, model.var), order = c(u = u, v = v, p = p, q = q), max.order = max.order, z = rep(0, times = length(series)), #這里的0，長度為444 h = rep(var(series), times = length(series)), #h是方差是1(因為此時的series是被標准化過的)，長度為444 x = series, scale = scale, init.rec = init.rec, h.start = h.start, llh.start = llh.start) ans}
其結果
這個函數的作用，主要是形成z、h之類的存儲容器
其中3
.params <- .garchInitParameters(formula.mean = formula.mean, formula.var = formula.var, delta = delta, skew = skew, shape = shape, cond.dist = cond.dist[1], include.mean = include.mean, include.delta = include.delta, include.skew = include.skew, include.shape = include.shape, leverage = leverage, algorithm = algorithm[1], control = con, trace = trace) .setfGarchEnv(.params = .params)
.garchInitParameters函數
function (formula.mean, formula.var, delta, skew, shape, cond.dist, include.mean, include.delta, include.skew, include.shape, leverage, algorithm, control, trace) { .DEBUG = FALSE .series <- .getfGarchEnv(".series") model.order = as.numeric(strsplit(strsplit(strsplit(as.character(formula.mean), "\\(")[[2]][2], "\\)")[[1]], ",")[[1]]) u = model.order[1] v = 0 if (length(model.order) == 2) v = model.order[2] model.order = as.numeric(strsplit(strsplit(strsplit(as.character(formula.var), "\\(")[[2]][2], "\\)")[[1]], ",")[[1]]) p = model.order[1] q = 0 if (length(model.order) == 2) q = model.order[2] model.var = .series$model[2] if (is.null(include.delta)) { if (model.var == "garch") { include.delta = FALSE } else { include.delta = TRUE } } if (is.null(leverage)) { if (model.var == "garch") { leverage = FALSE } else { leverage = TRUE } } if (cond.dist == "t") cond.dist = "std" skewed.dists = c("snorm", "sged", "sstd", "snig") if (is.null(include.skew)) { if (any(skewed.dists == cond.dist)) { include.skew = TRUE } else { include.skew = FALSE } } shaped.dists = c("ged", "sged", "std", "sstd", "snig") if (is.null(include.shape)) { if (any(shaped.dists == cond.dist)) { include.shape = TRUE } else { include.shape = FALSE } } Names = c("mu", if (u > 0) paste("ar", 1:u, sep = ""), if (v > 0) paste("ma", 1:v, sep = ""), "omega", if (p > 0) paste("alpha", 1:p, sep = ""), if (p > 0) paste("gamma", 1:p, sep = ""), if (q > 0) paste("beta", 1:q, sep = ""), "delta", "skew", "shape") fit.mean = arima(.series$x, order = c(u, 0, v), include.mean = include.mean)$coef alpha.start = 0.1 beta.start = 0.8 params = c( if (include.mean) fit.mean[length(fit.mean)] else 0, #最后一項就是intercept，所以就是mu啊！ if (u > 0) fit.mean[1:u], if (v > 0) fit.mean[(u + 1):(length(fit.mean) - as.integer(include.mean))], var(.series$x, na.rm = TRUE) * (1 - alpha.start - beta.start), if (p > 0) rep(alpha.start/p, times = p), if (p > 0) rep(0.1, times = p), if (q > 0) rep(beta.start/q, times = q), delta, skew, shape) names(params) = Names TINY = 1e-08 USKEW = 1/10 USHAPE = 1 if (cond.dist == "snig") USKEW = -0.99 U = c(-10 * abs(mean(.series$x)), if (u > 0) rep(-1 + TINY, times = u), if (v > 0) rep(-1 + TINY, times = v), 1e-06 * var(.series$x), if (p > 0) rep(0 + TINY, times = p), if (p > 0) rep(-1 + TINY, times = p), if (q > 0) rep(0 + TINY, times = q), 0, USKEW, USHAPE) names(U) = Names VSKEW = 10 VSHAPE = 10 if (cond.dist == "snig") VSKEW = 0.99 V = c(10 * abs(mean(.series$x)), if (u > 0) rep(1 - TINY, times = u), if (v > 0) rep(1 - TINY, times = v), 100 * var(.series$x), if (p > 0) rep(1 - TINY, times = p), if (p > 0) rep(1 - TINY, times = p), if (q > 0) rep(1 - TINY, times = q), 2, VSKEW, VSHAPE) names(V) = Names includes = c(include.mean, if (u > 0) rep(TRUE, times = u), if (v > 0) rep(TRUE, times = v), TRUE, if (p > 0) rep(TRUE, times = p), if (p > 0) rep(leverage, times = p), if (q > 0) rep(TRUE, times = q), include.delta, include.skew, include.shape) names(includes) = Names index = (1:length(params))[includes == TRUE] names(index) = names(params)[includes == TRUE] alpha <- beta <- NULL if (p > 0) alpha = params[substr(Names, 1, 5) == "alpha"] if (p > 0 & leverage) gamma = params[substr(Names, 1, 5) == "gamma"] if (p > 0 & !leverage) gamma = rep(0, times = p) if (q > 0) beta = params[substr(Names, 1, 4) == "beta"] if (.series$model[2] == "garch") { persistence = sum(alpha) + sum(beta) } else if (.series$model[2] == "aparch") { persistence = sum(beta) for (i in 1:p) persistence = persistence + alpha[i] * garchKappa(cond.dist, gamma[i], params["delta"], params["skew"], params["shape"]) } names(persistence) = "persistence" list(params = params, U = U, V = V, includes = includes, index = index, mu = params[1], delta = delta, skew = skew, shape = shape, cond.dist = cond.dist, leverage = leverage, persistence = persistence, control = control)}
這個函數的作用，主要是形成參數的初始值吧
fit.mean = arima(.series$x, order = c(u, 0, v), include.mean = include.mean)$coefparams = c( if (include.mean) fit.mean[length(fit.mean)] else 0
include.mean默認是T，fit.mean就相當於是對原始數據標准化后的序列進行arima建模，只不過ar和ma的階數都是0，只剩下截距項，作為mu的初始值，其中length(fit.mean)就是模型的截距項intercept。
此時的mu是 0.1128933
此時的mu是 0.1128933
事實上，如果是到底為止，那么就參數的值其實還是沒有被求出，為什么后面在
.garchFit
函數中卻直接開始賦值了呢？這是因為后面的.garchOptimizeLLH函數借助.get/.setfGarchEnv()函數，取得了.series和.prameter
其中4 (4.0.2.4).garchOptimizeLLH函數
.garchOptimizeLLH函數(刪去了一些用不到的代碼)
- function (hessian = hessian, robust.cvar, trace)
{ .series <- .getfGarchEnv(".series") .params <- .getfGarchEnv(".params") INDEX = .params$index algorithm = .params$control$algorithm[1] #algorithm是nlminb TOL1 = .params$control$tol1 TOL2 = .params$control$tol2 if (algorithm == "nlminb" | algorithm == "nlminb+nm") { fit <- .garchRnlminb(.params, .series, .garchLLH, trace) .params$llh = fit$llh .params$params[INDEX] = fit$par .setfGarchEnv(.params = .params) } .params$llh = fit$llh .params$params[INDEX] = fit$par .setfGarchEnv(.params = .params) if (hessian == "ropt") { fit$hessian <- -.garchRoptimhess(par = fit$par, .params = .params, .series = .series) titleHessian = "R-optimhess" }............... if (.params$control$xscale) { .series$x <- .series$x * .series$scale #將標准化后的值還原為原始的數值了 if (.params$include["mu"]) fit$coef["mu"] <- fit$par["mu"] <- .params$params["mu"] <- .params$params["mu"] * .series$scale if (.params$include["omega"]) fit$coef["omega"] <- fit$par["omega"] <- .params$params["omega"] <- .params$params["omega"] * .series$scale^(.params$params["delta"]) .setfGarchEnv(.params = .params) .setfGarchEnv(.series = .series) } if (.params$control$xscale) { if (.params$include["mu"]) { fit$hessian[, "mu"] <- fit$hessian[, "mu"]/.series$scale fit$hessian["mu", ] <- fit$hessian["mu", ]/.series$scale } if (.params$include["omega"]) { fit$hessian[, "omega"] <- fit$hessian[, "omega"]/.series$scale^(.params$params["delta"]) fit$hessian["omega", ] <- fit$hessian["omega", ]/.series$scale^(.params$params["delta"]) } } .llh <- fit$llh <- fit$value <- .garchLLH(fit$par, trace = FALSE, fGarchEnv = TRUE) .series <- .getfGarchEnv(".series") if (robust.cvar) fit$gradient <- -.garchRCDAGradient(par = fit$par, .params = .params, .series = .series) N = length(.series$x) NPAR = length(fit$par) fit$ics = c(AIC = c((2 * fit$value)/N + 2 * NPAR/N), BIC = (2 * fit$value)/N + NPAR * log(N)/N, SIC = (2 * fit$value)/N + log((N + 2 * NPAR)/N), HQIC = (2 * fit$value)/N + (2 * NPAR * log(log(N)))/N) names(fit$ics) <- c("AIC", "BIC", "SIC", "HQIC") fit}
其中4 (4.0.2.4.1)
.garchRnlminb函數
.garchRnlminb函數
function (.params, .series, .garchLLH, trace) { if (trace) cat("\nR coded nlminb Solver: \n\n") INDEX = .params$index parscale = rep(1, length = length(INDEX)) names(parscale) = names(.params$params[INDEX]) parscale["omega"] = var(.series$x)^(.params$delta/2) parscale["mu"] = abs(mean(.series$x)) #這個並非mu的初始值 TOL1 = .params$control$tol1 fit <- nlminb(start = .params$params[INDEX], objective = .garchLLH, lower = .params$U[INDEX], upper = .params$V[INDEX], scale = 1/parscale, control = list(eval.max = 2000, iter.max = 1500, rel.tol = 1e-14 * TOL1, x.tol = 1e-14 * TOL1, trace = as.integer(trace)), fGarchEnv = FALSE) fit$value <- fit.llh <- fit$objective names(fit$par) = names(.params$params[INDEX]) .garchLLH(fit$par, trace = FALSE, fGarchEnv = TRUE) fit$coef <- fit$par fit$llh <- fit$objective fit}
關於nlminb函數，請看筆記：
.params$params[INDEX]是初始值
objective = .garchLLH是被優化的函數
數據呢？怎么沒見傳入數據？
.garchRnlminb函數函數運行完，則會估計出參數：
fit$par
mu omega alpha1 beta1
0.08876900 0.05705989 0.08643831 0.85258554
這個值呢，是別標准化后的參數值
fit$coef["mu"] <- fit$par["mu"] <- .params$params["mu"] <- .params$params["mu"] * .series$scale
通過上述語句轉化為原始序列對應的參數即可！
Browse[4]> fit$coef
mu omega alpha1 beta1
0.0112656813 0.0009190163 0.0864383140 0.8525855370
這就是模型最終輸出的參數！
運行到.llh <- fit$llh <- fit$value <-
.garchLLH(fit$par, trace = FALSE,
fGarchEnv = TRUE)
這里跳進去
得到其中用到的參數.garchLLH函數
.garchLLHfunction (params, trace = TRUE, fGarchEnv = FALSE) { DEBUG = FALSE if (DEBUG) print("Entering Function .garchLLH") .series <- .getfGarchEnv(".series") .params <- .getfGarchEnv(".params") .garchDist <- .getfGarchEnv(".garchDist") #概率密度函數 .llh <- .getfGarchEnv(".llh") #對數似然函數 if (DEBUG) print(.params$control$llh) if (.params$control$llh == "internal") { if (DEBUG) print("internal") return(.aparchLLH.internal(params, trace = trace, fGarchEnv = fGarchEnv)) #結束 } else if (.params$control$llh == "filter") { if (DEBUG) print("filter") return(.aparchLLH.filter(params, trace = trace, fGarchEnv = fGarchEnv)) } else if (.params$control$llh == "testing") { if (DEBUG) print("testing") return(.aparchLLH.testing(params, trace = trace, fGarchEnv = fGarchEnv)) } else { stop("LLH is neither internal, testing, nor filter!") }}<environment: namespace:fGarch>
.aparchLLH.internal函數
其中我進去看了下
function (params, trace = TRUE, fGarchEnv = TRUE) { DEBUG = FALSE if (DEBUG) print("Entering Function .garchLLH.internal") .series <- .getfGarchEnv(".series") .params <- .getfGarchEnv(".params") .garchDist <- .getfGarchEnv(".garchDist") .llh <- .getfGarchEnv(".llh") if (DEBUG) print(.params$control$llh) if (.params$control$llh == "internal") { INDEX <- .params$index MDIST <- c(norm = 10, QMLE = 10, snorm = 11, std = 20, sstd = 21, ged = 30, sged = 31)[.params$cond.dist] if (.params$control$fscale) NORM <- length(.series$x) else NORM = 1 REC <- 1 if (.series$init.rec == "uev") REC <- 2 MYPAR <- c(REC = REC, LEV = as.integer(.params$leverage), MEAN = as.integer(.params$includes["mu"]), DELTA = as.integer(.params$includes["delta"]), SKEW = as.integer(.params$includes["skew"]), SHAPE = as.integer(.params$includes["shape"]), ORDER = .series$order, NORM = as.integer(NORM)) MAX <- max(.series$order) NF <- length(INDEX) N <- length(.series$x) DPARM <- c(.params$delta, .params$skew, .params$shape) fit <- .Fortran("garchllh", N = as.integer(N), Y = as.double(.series$x), Z = as.double(.series$z), H = as.double(.series$h), NF = as.integer(NF), X = as.double(params), DPARM = as.double(DPARM), MDIST = as.integer(MDIST), MYPAR = as.integer(MYPAR), F = as.double(0), PACKAGE = "fGarch") llh <- fit[[10]] #此時llh就是-312 if (is.na(llh)) llh = .llh + 0.1 * (abs(.llh)) if (!is.finite(llh)) llh = .llh + 0.1 * (abs(.llh)) .setfGarchEnv(.llh = llh) if (fGarchEnv) { .series$h <- fit[[4]] .series$z <- fit[[3]] #直到這里，z才從444個0被替換為殘差 .setfGarchEnv(.series = .series) } } else { stop("LLH is not internal!") } c(LogLikelihood = llh)}
.Fortran是調用Fortran代碼的，真是我擦！
Browse[5]> .garchDistfunction (z, hh, skew, shape) { dnorm(x = z/hh, mean = 0, sd = 1)/hh}<environment: namespace:fGarch>
寫在2017年1月28日00:54:55
農歷新年凌晨，新的一年要更加努力~