使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码

时间：2021-09-02 23:03:00

一、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.

二、查看数据----------------------------------------------
使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码
双击da(相当于View(da))，可以看到数据是这样的
使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码

经过对数转化后的收益率，是我们建模用的原始数据。
    

    2. 

  2. intc=log(da$intc+1)
    3. 

  3. intc
    4. 

  4.   [1]  0.0099998346 -0.1500127529  0.0670640792
    5. 

  5.   [4]  0.0829486352 -0.1103484905  0.1251628488
    6. 

  6.   [7]  0.4855078158  0.1112255825  0.2109235908
    7. 

  7. ..........
    8. 


    

  1. library(fGarch)
    2. 

  2. da=read.table("m-intcsp7309.txt",header=T)
    3. 

  3. intc=log(da$intc+1)
    4. 

  4. m4=garchFit(~1+garch(1,1),data=intc)
    5. 

  5. summary(m4)
    6. 

  6. Title:
    7. 

  7.  GARCH Modelling 
    8. 

    9. 

  9. Call:
    10. 

  10.  garchFit(formula = ~1 + garch(1, 1), data = intc) 
    11. 

    12. 

  12. Mean and Variance Equation:
    13. 

  13.  data ~ 1 + garch(1, 1)
    14. 

  14. <environment: 0x000000001d475658>
    15. 

  15.  [data = intc]
    16. 

    17. 

  17. Conditional Distribution:
    18. 

  18.  norm 
    19. 

    20. 

  20. Coefficient(s):
    21. 

  21.         mu       omega      alpha1       beta1  
    22. 

  22. 0.01126568  0.00091902  0.08643831  0.85258554  
    23. 

    24. 

  24. Std. Errors:
    25. 

  25.  based on Hessian 
    26. 

    27. 

  27. Error Analysis:
    28. 

  28.         Estimate  Std. Error  t value Pr(>|t|)    
    29. 

  29. mu     0.0112657   0.0053931    2.089  0.03672 *  
    30. 

  30. omega  0.0009190   0.0003888    2.364  0.01808 *  
    31. 

  31. alpha1 0.0864383   0.0265439    3.256  0.00113 ** 
    32. 

  32. beta1  0.8525855   0.0394322   21.622  < 2e-16 ***
    33. 

  33. ---
    34. 

  34. Signif. codes:  
    35. 

  35. 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    36. 

    37. 

  37. Log Likelihood:
    38. 

  38.  312.3307    normalized:  0.7034475 
    39. 

    40. 

  40. Description:
    41. 

  41.  Fri Jan 27 21:49:06 2017 by user: xuan 
    42. 

    43. 

    44. 

  44. Standardised Residuals Tests:
    45. 

  45.                                 Statistic p-Value     
    46. 

  46.  Jarque-Bera Test   R    Chi^2  174.904   0           
    47. 

  47.  Shapiro-Wilk Test  R    W      0.9709615 1.030282e-07
    48. 

  48.  Ljung-Box Test     R    Q(10)  8.016844  0.6271916   
    49. 

  49.  Ljung-Box Test     R    Q(15)  15.5006   0.4159946   
    50. 

  50.  Ljung-Box Test     R    Q(20)  16.41549  0.6905368   
    51. 

  51.  Ljung-Box Test     R^2  Q(10)  0.8746345 0.9999072   
    52. 

  52.  Ljung-Box Test     R^2  Q(15)  11.35935  0.7267295   
    53. 

  53.  Ljung-Box Test     R^2  Q(20)  12.55994  0.8954573   
    54. 

  54.  LM Arch Test       R    TR^2   10.51401  0.5709617   
    55. 

    56. 

  56. Information Criterion Statistics:
    57. 

  57.       AIC       BIC       SIC      HQIC 
    58. 

  58. -1.388877 -1.351978 -1.389037 -1.374326 
    59.


三、主要的函数结构----------------------------------------------
garchFit函数
         match.xxx函数....
         各种判断可以忽略
         .garchArgsParser
         .garchFit函数    
         ans最终的返回结果
         
         .garchFit函数
                       .garchInitSeries函数
                       .garchInitParameters函数
                       .garchOptimizeLLH函数
                                    .garchRnlminb函数
                                              nlminb函数
                                    .garchLLH函数
                                              .Fortran("garchllh")


四、调试----------------------------------------------

(4.0)garchFit函数
    

  1. function (formula = ~garch(1, 1), data = dem2gbp, init.rec = c("mci", 
    2. 

  2.   "uev"), delta = 2, skew = 1, shape = 4, cond.dist = c("norm", 
    3. 

  3.   "snorm", "ged", "sged", "std", "sstd", "snig", "QMLE"), 
    4. 

  4.   include.mean = TRUE, include.delta = NULL, include.skew = NULL, 
    5. 

  5.   include.shape = NULL, leverage = NULL, trace = TRUE, algorithm = c("nlminb", 
    6. 

  6.     "lbfgsb", "nlminb+nm", "lbfgsb+nm"), hessian = c("ropt", 
    7. 

  7.     "rcd"), control = list(), title = NULL, description = NULL, 
    8. 

  8.   ...) 
    9. 

  9. {
    10. 

  10.   DEBUG = FALSE
    11. 

  11.   init.rec = match.arg(init.rec)
    12. 

  12.   cond.dist = match.arg(cond.dist)
    13. 

  13.   hessian = match.arg(hessian)
    14. 

  14.   algorithm = match.arg(algorithm)
    15. 

  15.   CALL = match.call()
    16. 

  16.   Name = capture.output(substitute(data))
    17. 

  17.   if (is.character(data)) {
    18. 

  18.     eval(parse(text = paste("data(", data, ")")))
    19. 

  19.     data = eval(parse(text = data))
    20. 

  20.   }
    21. 

  21.   data <- as.data.frame(data)
    22. 

  22.   if (isUnivariate(data)) {
    23. 

  23.     colnames(data) <- "data"
    24. 

  24.   }
    25. 

  25.   else {
    26. 

  26.     uniqueNames = unique(sort(colnames(data)))
    27. 

  27.     if (is.null(colnames(data))) {
    28. 

  28.       stop("Column names of data are missing.")
    29. 

  29.     }
    30. 

  30.     if (length(colnames(data)) != length(uniqueNames)) {
    31. 

  31.       stop("Column names of data are not unique.")
    32. 

  32.     }
    33. 

  33.   }
    34. 

  34.   if (length(formula) == 3 && isUnivariate(data)) 
    35. 

  35.     formula[2] <- NULL
    36. 

  36.   if (length(formula) == 2) {
    37. 

  37.     if (isUnivariate(data)) {
    38. 

  38.       formula = as.formula(paste("data", paste(formula, 
    39. 

  39.         collapse = " ")))
    40. 

  40.     }
    41. 

  41.     else {
    42. 

  42.       stop("Multivariate data inputs require lhs for the formula.")
    43. 

  43.     }
    44. 

  44.   }
    45. 

  45.   robust.cvar <- (cond.dist == "QMLE")
    46. 

  46.   args = .garchArgsParser(formula = formula, data = data, 
    47. 

  47.     trace = FALSE)
    48. 

  48.   if (DEBUG) 
    49. 

  49.     print(list(formula.mean = args$formula.mean, formula.var = args$formula.var, 
    50. 

  50.       series = args$series, init.rec = init.rec, delta = delta, 
    51. 

  51.       skew = skew, shape = shape, cond.dist = cond.dist, 
    52. 

  52.       include.mean = include.mean, include.delta = include.delta, 
    53. 

  53.       include.skew = include.skew, include.shape = include.shape, 
    54. 

  54.       leverage = leverage, trace = trace, algorithm = algorithm, 
    55. 

  55.       hessian = hessian, robust.cvar = robust.cvar, control = control, 
    56. 

  56.       title = title, description = description))
    57. 

  57.   ans = .garchFit(formula.mean = args$formula.mean, formula.var = args$formula.var, 
    58. 

  58.     series = args$series, init.rec, delta, skew, shape, 
    59. 

  59.     cond.dist, include.mean, include.delta, include.skew, 
    60. 

  60.     include.shape, leverage, trace, algorithm, hessian, 
    61. 

  61.     robust.cvar, control, title, description, ...)
    62. 

  62.   ans@call = CALL
    63. 

  63.   attr(formula, "data") <- paste("data = ", Name, sep = "")
    64. 

  64.   ans@formula = formula
    65. 

  65.   ans
    66. 

  66. }
    67.

以下四句都是对传入的参数是否math(合法、规范)做检验的
    

  1.   init.rec = match.arg(init.rec)
    2. 

  2.   cond.dist = match.arg(cond.dist)
    3. 

  3.   hessian = match.arg(hessian)
    4. 

  4.   algorithm = match.arg(algorithm)
    5. 

  5.   CALL = match.call()
    6.
我们可以按下使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码进入函数查看
截取match.arg函数中的几句作说明：
    

  1. formal.args <- formals(sys.function(sys.parent()))
    2. 

  2. .......................................
    3. 

  3. if (identical(arg, choices)) 
    4. 

  4.       return(arg[1L])
    5.

formal argument形参

如果传入和默认的一样，说明没有传入，就返回第一个作为默认值
可以点击使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码跳出函数
 

(4.0.1).garchArgsParser函数
  1. args = .garchArgsParser(formula = formula, data = data, trace = FALSE)
    

  1. function (formula, data, trace = FALSE) 
    2. 

  2. {
    3. 

  3.     allVars = unique(sort(all.vars(formula)))
    4. 

  4.     allVarsTest = mean(allVars %in% colnames(data))
    5. 

  5.     if (allVarsTest != 1) {
    6. 

  6.         print(allVars)
    7. 

  7.         print(colnames(data))
    8. 

  8.         stop("Formula and data units do not match.")
    9. 

  9.     }
    10. 

  10.     formula.lhs = as.character(formula)[2]
    11. 

  11.     mf = match.call(expand.dots = FALSE)
    12. 

  12.     if (trace) {
    13. 

  13.         cat("\nMatched Function Call:\n ")
    14. 

  14.         print(mf)
    15. 

  15.     }
    16. 

  16.     m = match(c("formula", "data"), names(mf), 0)
    17. 

  17.     mf = mf[c(1, m)]
    18. 

  18.     mf[[1]] = as.name(".garchModelSeries")
    19. 

  19.     mf$fake = FALSE
    20. 

  20.     mf$lhs = TRUE
    21. 

  21.     if (trace) {
    22. 

  22.         cat("\nModelSeries Call:\n ")
    23. 

  23.         print(mf)
    24. 

  24.     }
    25. 

  25.     x = eval(mf, parent.frame())
    26. 

  26.     if (trace) 
    27. 

  27.         print(x)
    28. 

  28.     x = as.vector(x[, 1])
    29. 

  29.     names(x) = rownames(data)
    30. 

  30.     if (trace) 
    31. 

  31.         print(x)
    32. 

  32.     allLabels = attr(terms(formula), "term.labels")
    33. 

  33.     if (trace) {
    34. 

  34.         cat("\nAll Term Labels:\n ")
    35. 

  35.         print(allLabels)
    36. 

  36.     }
    37. 

  37.     if (length(allLabels) == 2) {
    38. 

  38.         formula.mean = as.formula(paste("~", allLabels[1]))
    39. 

  39.         formula.var = as.formula(paste("~", allLabels[2]))
    40. 

  40.     }
    41. 

  41.     else if (length(allLabels) == 1) {
    42. 

  42.         formula.mean = as.formula("~ arma(0, 0)")
    43. 

  43.         formula.var = as.formula(paste("~", allLabels[1]))
    44. 

  44.     }
    45. 

  45.     if (trace) {
    46. 

  46.         cat("\nMean Formula:\n ")
    47. 

  47.         print(formula.mean)
    48. 

  48.         cat("\nVariance Formula:\n ")
    49. 

  49.         print(formula.var)
    50. 

  50.     }
    51. 

  51.     ans <- list(formula.mean = formula.mean, formula.var = formula.var, 
    52. 

  52.         formula.lhs = formula.lhs, series = x)
    53. 

  53.     ans
    54. 

  54. }
    55.
提下这个函数里面的一些信息:
data是数据
%in%运算符就是math的作用 ?'%in%'
x就是原始数据
(需要特别说明的是,x是从data里取出来后as.vector的，被转化为向量后，会四舍五入显示了，所以如果直接在environment pane中查看有点不橡原始数据，比如intc的第一个值是0.0099998346，而这里x的值显示是0.01)
函数运行到最后一行的时候，environment pane中的部分信息。
使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码(方程的基本形式已经被解析出来了)
使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码

(4.0.2).garchFit函数
    

  1. ans = .garchFit(formula.mean = args$formula.mean, formula.var = args$formula.var, 
    2. 

  2.     series = args$series, init.rec, delta, skew, shape, 
    3. 

  3.     cond.dist, include.mean, include.delta, include.skew, 
    4. 

  4.     include.shape, leverage, trace, algorithm, hessian, 
    5. 

  5.     robust.cvar, control, title, description, ...)
    6.
我们为什么分析了.garchArgsParser这个解析函数，就是因为.garchFit函数用到的许多参数都是从args中取出的。
另外提到一下，带.开头的函数和变量都是内部的，不会直接显示出来的。
如果是变量名.xxx，这样的形式，那么需要我们自己在控制台输入.xxx查看，不能在envir pane中查看。
如果要查看内部函数，则
使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码
单击可以切换到对应的源代码和变量
使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码可以选择环境，一般情况是在global environment
.garchFit函数的源码
    

  1. function (formula.mean = ~arma(0, 0), formula.var = ~garch(1, 
    2. 

  2.   1), series, init.rec = c("mci", "uev"), delta = 2, skew = 1, 
    3. 

  3.   shape = 4, cond.dist = c("norm", "snorm", "ged", "sged", 
    4. 

  4.     "std", "sstd", "QMLE"), include.mean = TRUE, include.delta = NULL, 
    5. 

  5.   include.skew = NULL, include.shape = NULL, leverage = NULL, 
    6. 

  6.   trace = TRUE, algorithm = c("sqp", "nlminb", "lbfgsb", "nlminb+nm", 
    7. 

  7.     "lbfgsb+nm"), hessian = c("ropt", "rcd"), robust.cvar, 
    8. 

  8.   control = list(), title = NULL, description = NULL, ...) 
    9. 

  9. {
    10. 

  10.   DEBUG <- FALSE
    11. 

  11.   if (DEBUG) 
    12. 

  12.     print("Formula Specification ...")
    13. 

  13.   fcheck = rev(all.names(formula.mean))[1]
    14. 

  14.   if (fcheck == "ma") {
    15. 

  15.     stop("Use full formula: arma(0,q) for ma(q)")
    16. 

  16.   }
    17. 

  17.   else if (fcheck == "ar") {
    18. 

  18.     stop("Use full formula expression: arma(p,0) for ar(p)")
    19. 

  19.   }
    20. 

  20.   if (DEBUG) 
    21. 

  21.     print("Recursion Initialization ...")
    22. 

  22.   if (init.rec[1] != "mci" & algorithm[1] != "sqp") {
    23. 

  23.     stop("Algorithm only supported for mci Recursion")
    24. 

  24.   }
    25. 

  25.   .StartFit <- Sys.time()
    26. 

  26.   if (DEBUG) 
    27. 

  27.     print("Generate Control List ...")
    28. 

  28.   con <- .garchOptimizerControl(algorithm, cond.dist)
    29. 

  29.   con[(namc <- names(control))] <- control
    30. 

  30.   if (DEBUG) 
    31. 

  31.     print("Initialize Time Series ...")
    32. 

  32.   data <- series
    33. 

  33.   scale <- if (con$xscale) 
    34. 

  34.     sd(series)
    35. 

  35.   else 1
    36. 

  36.   series <- series/scale
    37. 

  37.   .series <- .garchInitSeries(formula.mean = formula.mean, 
    38. 

  38.     formula.var = formula.var, cond.dist = cond.dist[1], 
    39. 

  39.     series = series, scale = scale, init.rec = init.rec[1], 
    40. 

  40.     h.start = NULL, llh.start = NULL, trace = trace)
    41. 

  41.   .setfGarchEnv(.series = .series)
    42. 

  42.   if (DEBUG) 
    43. 

  43.     print("Initialize Model Parameters ...")
    44. 

  44.   .params <- .garchInitParameters(formula.mean = formula.mean, 
    45. 

  45.     formula.var = formula.var, delta = delta, skew = skew, 
    46. 

  46.     shape = shape, cond.dist = cond.dist[1], include.mean = include.mean, 
    47. 

  47.     include.delta = include.delta, include.skew = include.skew, 
    48. 

  48.     include.shape = include.shape, leverage = leverage, 
    49. 

  49.     algorithm = algorithm[1], control = con, trace = trace)
    50. 

  50.   .setfGarchEnv(.params = .params)
    51. 

  51.   if (DEBUG) 
    52. 

  52.     print("Select Conditional Distribution ...")
    53. 

  53.   .setfGarchEnv(.garchDist = .garchSetCondDist(cond.dist[1]))
    54. 

  54.   if (DEBUG) 
    55. 

  55.     print("Estimate Model Parameters ...")
    56. 

  56.   .setfGarchEnv(.llh = 1e+99)
    57. 

  57.   .llh <- .getfGarchEnv(".llh")
    58. 

  58.   fit = .garchOptimizeLLH(hessian, robust.cvar, trace)
    59. 

  59.   if (DEBUG) 
    60. 

  60.     print("Add to fit ...")
    61. 

  61.   .series <- .getfGarchEnv(".series")
    62. 

  62.   .params <- .getfGarchEnv(".params")
    63. 

  63.   names(.series$h) <- NULL
    64. 

  64.   fit$series = .series
    65. 

  65.   fit$params = .params
    66. 

  66.   if (DEBUG) 
    67. 

  67.     print("Retrieve Residuals and Fitted Values ...")
    68. 

  68.   residuals = .series$z
    69. 

  69.   fitted.values = .series$x - residuals
    70. 

  70.   h.t = .series$h
    71. 

  71.   if (.params$includes["delta"]) 
    72. 

  72.     deltainv = 1/fit$par["delta"]
    73. 

  73.   else deltainv = 1/fit$params$delta
    74. 

  74.   sigma.t = (.series$h)^deltainv
    75. 

  75.   if (DEBUG) 
    76. 

  76.     print("Standard Errors and t-Values ...")
    77. 

  77.   fit$cvar <- if (robust.cvar) 
    78. 

  78.     (solve(fit$hessian) %*% (t(fit$gradient) %*% fit$gradient) %*% 
    79. 

  79.       solve(fit$hessian))
    80. 

  80.   else -solve(fit$hessian)
    81. 

  81.   fit$se.coef = sqrt(diag(fit$cvar))
    82. 

  82.   fit$tval = fit$coef/fit$se.coef
    83. 

  83.   fit$matcoef = cbind(fit$coef, fit$se.coef, fit$tval, 2 * 
    84. 

  84.     (1 - pnorm(abs(fit$tval))))
    85. 

  85.   dimnames(fit$matcoef) = list(names(fit$tval), c(" Estimate", 
    86. 

  86.     " Std. Error", " t value", "Pr(>|t|)"))
    87. 

  87.   if (DEBUG) 
    88. 

  88.     print("Add Title and Description ...")
    89. 

  89.   if (is.null(title)) 
    90. 

  90.     title = "GARCH Modelling"
    91. 

  91.   if (is.null(description)) 
    92. 

  92.     description = description()
    93. 

  93.   Time = Sys.time() - .StartFit
    94. 

  94.   if (trace) {
    95. 

  95.     cat("\nTime to Estimate Parameters:\n ")
    96. 

  96.     print(Time)
    97. 

  97.   }
    98. 

  98.   new("fGARCH", call = as.call(match.call()), formula = as.formula(paste("~", 
    99. 

  99.     formula.mean, "+", formula.var)), method = "Max Log-Likelihood Estimation", 
    100. 

  100.     data = data, fit = fit, residuals = residuals, fitted = fitted.values, 
    101. 

  101.     h.t = h.t, sigma.t = as.vector(sigma.t), title = as.character(title), 
    102. 

  102.     description = as.character(description))
    103. 

  103. }
    104.

在正式分析函数的源码之前
先看一下模型变量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
    

  1. scale <- if (con$xscale) 
    2. 

  2.     sd(series)
    3. 

  3. series <- series/scale
    4.
是对数据进行标准化
（注意，标准化之后，数据和原始数据不同啦，在后面运算时，不要忘了~）
我们可以在调试界面的控制台直接输入命令，来查看变量在当前环境和一定语句代运行完后的值
使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码
这个所谓的标准化，并没有减去均值的，要注意！

其中2 
    

  1. .series <- .garchInitSeries(formula.mean = formula.mean, 
    2. 

  2.     formula.var = formula.var, cond.dist = cond.dist[1], 
    3. 

  3.     series = series, scale = scale, init.rec = init.rec[1], 
    4. 

  4.     h.start = NULL, llh.start = NULL, trace = trace)
    5. 

  5. .setfGarchEnv(.series = .series)
    6.
注意这种.setfGarchEnv的做法，就是对环境变量中的变量赋值
还有.getfGarchEnv，可以提高自己的编程能力。
.garchInitSeries函数
    

  1. function (formula.mean, formula.var, cond.dist, series, scale, 
    2. 

  2.   init.rec, h.start, llh.start, trace) 
    3. 

  3. {
    4. 

  4.   mm = length(formula.mean)
    5. 

  5.   if (mm != 2) 
    6. 

  6.     stop("Mean Formula misspecified")
    7. 

  7.   end = regexpr("\\(", as.character(formula.mean[mm])) - 1
    8. 

  8.   model.mean = substr(as.character(formula.mean[mm]), 1, end)
    9. 

  9.   if (!any(c("ar", "ma", "arma") == model.mean)) 
    10. 

  10.     stop("formula.mean must be one of: ar, ma, arma")
    11. 

  11.   mv = length(formula.var)
    12. 

  12.   if (mv != 2) 
    13. 

  13.     stop("Variance Formula misspecified")
    14. 

  14.   end = regexpr("\\(", as.character(formula.var[mv])) - 1
    15. 

  15.   model.var = substr(as.character(formula.var[mv]), 1, end)
    16. 

  16.   if (!any(c("garch", "aparch") == model.var)) 
    17. 

  17.     stop("formula.var must be one of: garch, aparch")
    18. 

  18.   model.order = as.numeric(strsplit(strsplit(strsplit(as.character(formula.mean), 
    19. 

  19.     "\\(")[[2]][2], "\\)")[[1]], ",")[[1]])
    20. 

  20.   u = model.order[1]
    21. 

  21.   v = 0
    22. 

  22.   if (length(model.order) == 2) 
    23. 

  23.     v = model.order[2]
    24. 

  24.   maxuv = max(u, v)
    25. 

  25.   if (u < 0 | v < 0) 
    26. 

  26.     stop("*** ARMA orders must be positive.")
    27. 

  27.   model.order = as.numeric(strsplit(strsplit(strsplit(as.character(formula.var), 
    28. 

  28.     "\\(")[[2]][2], "\\)")[[1]], ",")[[1]])
    29. 

  29.   p = model.order[1]
    30. 

  30.   q = 0
    31. 

  31.   if (length(model.order) == 2) 
    32. 

  32.     q = model.order[2]
    33. 

  33.   if (p + q == 0) 
    34. 

  34.     stop("Misspecified GARCH Model: Both Orders are zero!")
    35. 

  35.   maxpq = max(p, q)
    36. 

  36.   if (p < 0 | q < 0) 
    37. 

  37.     stop("*** GARCH orders must be positive.")
    38. 

  38.   max.order = max(maxuv, maxpq)
    39. 

  39.   if (is.null(h.start)) 
    40. 

  40.     h.start = max.order + 1
    41. 

  41.   if (is.null(llh.start)) 
    42. 

  42.     llh.start = 1
    43. 

  43.   if (init.rec != "mci" & model.var != "garch") {
    44. 

  44.     stop("Algorithm only supported for mci Recursion")
    45. 

  45.   }
    46. 

  46.  
    47. 

  47.   ans = list(model = c(model.mean, model.var), order = c(u = u, 
    48. 

  48.     v = v, p = p, q = q), max.order = max.order,
    49. 

  49.     z = rep(0, times = length(series)),   #这里的0，长度为444
    50. 

  50.     h = rep(var(series), times = length(series)),  #h是方差是1(因为此时的series是被标准化过的)，长度为444
    51. 

  51.     x = series, scale = scale, init.rec = init.rec, h.start = h.start, 
    52. 

  52.     llh.start = llh.start)
    53. 

  53.   ans
    54. 

  54. }
    55.

其结果
使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码
 使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码
 使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码
 使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码
这个函数的作用，主要是形成z、h之类的存储容器


其中3
    

  1. .params <- .garchInitParameters(formula.mean = formula.mean, 
    2. 

  2.     formula.var = formula.var, delta = delta, skew = skew, 
    3. 

  3.     shape = shape, cond.dist = cond.dist[1], include.mean = include.mean, 
    4. 

  4.     include.delta = include.delta, include.skew = include.skew, 
    5. 

  5.     include.shape = include.shape, leverage = leverage, 
    6. 

  6.     algorithm = algorithm[1], control = con, trace = trace)
    7. 

  7.   .setfGarchEnv(.params = .params)
    8.

.garchInitParameters函数 
    

  1. function (formula.mean, formula.var, delta, skew, shape, cond.dist, 
    2. 

  2.   include.mean, include.delta, include.skew, include.shape, 
    3. 

  3.   leverage, algorithm, control, trace) 
    4. 

  4. {
    5. 

  5.   .DEBUG = FALSE
    6. 

  6.   .series <- .getfGarchEnv(".series")
    7. 

  7.   model.order = as.numeric(strsplit(strsplit(strsplit(as.character(formula.mean), 
    8. 

  8.     "\\(")[[2]][2], "\\)")[[1]], ",")[[1]])
    9. 

  9.   u = model.order[1]
    10. 

  10.   v = 0
    11. 

  11.   if (length(model.order) == 2) 
    12. 

  12.     v = model.order[2]
    13. 

  13.   model.order = as.numeric(strsplit(strsplit(strsplit(as.character(formula.var), 
    14. 

  14.     "\\(")[[2]][2], "\\)")[[1]], ",")[[1]])
    15. 

  15.   p = model.order[1]
    16. 

  16.   q = 0
    17. 

  17.   if (length(model.order) == 2) 
    18. 

  18.     q = model.order[2]
    19. 

  19.   model.var = .series$model[2]
    20. 

  20.   if (is.null(include.delta)) {
    21. 

  21.     if (model.var == "garch") {
    22. 

  22.       include.delta = FALSE
    23. 

  23.     }
    24. 

  24.     else {
    25. 

  25.       include.delta = TRUE
    26. 

  26.     }
    27. 

  27.   }
    28. 

  28.   if (is.null(leverage)) {
    29. 

  29.     if (model.var == "garch") {
    30. 

  30.       leverage = FALSE
    31. 

  31.     }
    32. 

  32.     else {
    33. 

  33.       leverage = TRUE
    34. 

  34.     }
    35. 

  35.   }
    36. 

  36.   if (cond.dist == "t") 
    37. 

  37.     cond.dist = "std"
    38. 

  38.   skewed.dists = c("snorm", "sged", "sstd", "snig")
    39. 

  39.   if (is.null(include.skew)) {
    40. 

  40.     if (any(skewed.dists == cond.dist)) {
    41. 

  41.       include.skew = TRUE
    42. 

  42.     }
    43. 

  43.     else {
    44. 

  44.       include.skew = FALSE
    45. 

  45.     }
    46. 

  46.   }
    47. 

  47.   shaped.dists = c("ged", "sged", "std", "sstd", "snig")
    48. 

  48.   if (is.null(include.shape)) {
    49. 

  49.     if (any(shaped.dists == cond.dist)) {
    50. 

  50.       include.shape = TRUE
    51. 

  51.     }
    52. 

  52.     else {
    53. 

  53.       include.shape = FALSE
    54. 

  54.     }
    55. 

  55.   }
    56. 

  56.   Names = c("mu", 
    57. 

  57.                      if (u > 0) paste("ar", 1:u, sep = ""), 
    58. 

  58.                      if (v > 0) paste("ma", 1:v, sep = ""), "omega", 
    59. 

  59.                      if (p > 0) paste("alpha", 1:p, sep = ""), 
    60. 

  60.                      if (p > 0) paste("gamma", 1:p, sep = ""), 
    61. 

  61.                      if (q > 0) paste("beta", 1:q, sep = ""), 
    62. 

  62.                      "delta", "skew", "shape")
    63. 

  63.   fit.mean = arima(.series$x, order = c(u, 0, v), include.mean = include.mean)$coef
    64. 

  64.   alpha.start = 0.1
    65. 

  65.   beta.start = 0.8
    66. 

  66.   params = c(
    67. 

  67.     if (include.mean) fit.mean[length(fit.mean)] else 0,      #最后一项就是intercept，所以就是mu啊！
    68. 

  68.     if (u > 0) fit.mean[1:u], 
    69. 

  69.     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), 
    70. 

  70.     if (p > 0) rep(alpha.start/p, times = p), 
    71. 

  71.     if (p > 0) rep(0.1, times = p), 
    72. 

  72.     if (q > 0) rep(beta.start/q, times = q), 
    73. 

  73.     delta, skew, shape)
    74. 

  74.   names(params) = Names
    75. 

  75.   TINY = 1e-08
    76. 

  76.   USKEW = 1/10
    77. 

  77.   USHAPE = 1
    78. 

  78.   if (cond.dist == "snig") 
    79. 

  79.     USKEW = -0.99
    80. 

  80.   U = c(-10 * abs(mean(.series$x)), if (u > 0) rep(-1 + TINY, 
    81. 

  81.     times = u), if (v > 0) rep(-1 + TINY, times = v), 1e-06 * 
    82. 

  82.     var(.series$x), if (p > 0) rep(0 + TINY, times = p), 
    83. 

  83.     if (p > 0) rep(-1 + TINY, times = p), if (q > 0) rep(0 + 
    84. 

  84.       TINY, times = q), 0, USKEW, USHAPE)
    85. 

  85.   names(U) = Names
    86. 

  86.   VSKEW = 10
    87. 

  87.   VSHAPE = 10
    88. 

  88.   if (cond.dist == "snig") 
    89. 

  89.     VSKEW = 0.99
    90. 

  90.   V = c(10 * abs(mean(.series$x)), if (u > 0) rep(1 - TINY, 
    91. 

  91.     times = u), if (v > 0) rep(1 - TINY, times = v), 100 * 
    92. 

  92.     var(.series$x), if (p > 0) rep(1 - TINY, times = p), 
    93. 

  93.     if (p > 0) rep(1 - TINY, times = p), if (q > 0) rep(1 - 
    94. 

  94.       TINY, times = q), 2, VSKEW, VSHAPE)
    95. 

  95.   names(V) = Names
    96. 

  96.   includes = c(include.mean, if (u > 0) rep(TRUE, times = u), 
    97. 

  97.     if (v > 0) rep(TRUE, times = v), TRUE, if (p > 0) rep(TRUE, 
    98. 

  98.       times = p), if (p > 0) rep(leverage, times = p), 
    99. 

  99.     if (q > 0) rep(TRUE, times = q), include.delta, include.skew, 
    100. 

  100.     include.shape)
    101. 

  101.   names(includes) = Names
    102. 

  102.   index = (1:length(params))[includes == TRUE]
    103. 

  103.   names(index) = names(params)[includes == TRUE]
    104. 

  104.   alpha <- beta <- NULL
    105. 

  105.   if (p > 0) 
    106. 

  106.     alpha = params[substr(Names, 1, 5) == "alpha"]
    107. 

  107.   if (p > 0 & leverage) 
    108. 

  108.     gamma = params[substr(Names, 1, 5) == "gamma"]
    109. 

  109.   if (p > 0 & !leverage) 
    110. 

  110.     gamma = rep(0, times = p)
    111. 

  111.   if (q > 0) 
    112. 

  112.     beta = params[substr(Names, 1, 4) == "beta"]
    113. 

  113.   if (.series$model[2] == "garch") {
    114. 

  114.     persistence = sum(alpha) + sum(beta)
    115. 

  115.   }
    116. 

  116.   else if (.series$model[2] == "aparch") {
    117. 

  117.     persistence = sum(beta)
    118. 

  118.     for (i in 1:p) 
    119. 

  119.           persistence = persistence + alpha[i] * garchKappa(cond.dist, gamma[i], params["delta"], params["skew"], params["shape"])
    120. 

  120.   }
    121. 

  121.   names(persistence) = "persistence"
    122. 

  122.   list(params = params, U = U, V = V, includes = includes, 
    123. 

  123.     index = index, mu = params[1], delta = delta, skew = skew, 
    124. 

  124.     shape = shape, cond.dist = cond.dist, leverage = leverage, 
    125. 

  125.     persistence = persistence, control = control)
    126. 

  126. }
    127.
这个函数的作用，主要是形成参数的初始值吧
    

  1. fit.mean = arima(.series$x, order = c(u, 0, v), include.mean = include.mean)$coef
    2. 

  2. params = c(
    3. 

  3.     if (include.mean) fit.mean[length(fit.mean)] else 0
    4.
include.mean默认是T，fit.mean就相当于是对原始数据标准化后的序列进行arima建模，只不过ar和ma的阶数都是0，只剩下截距项，作为mu的初始值，其中length(fit.mean)就是模型的截距项intercept。
此时的mu是0.1128933


事实上，如果是到底为止，那么就参数的值其实还是没有被求出，为什么后面在.garchFit函数中却直接开始赋值了呢？这是因为后面的.garchOptimizeLLH函数借助.get/.setfGarchEnv()函数，取得了.series和.prameter


其中4 (4.0.2.4).garchOptimizeLLH函数
.garchOptimizeLLH函数(删去了一些用不到的代码)
    

  1. function (hessian = hessian, robust.cvar, trace)
    2. 

  2. {
    3. 

  3.   .series <- .getfGarchEnv(".series")
    4. 

  4.   .params <- .getfGarchEnv(".params")
    5. 

  5.   INDEX = .params$index
    6. 

  6.   algorithm = .params$control$algorithm[1]   #algorithm是nlminb
    7. 

  7.   TOL1 = .params$control$tol1
    8. 

  8.   TOL2 = .params$control$tol2
    9. 

  9.   if (algorithm == "nlminb" | algorithm == "nlminb+nm") {
    10. 

  10.     fit <- .garchRnlminb(.params, .series, .garchLLH, trace)
    11. 

  11.     .params$llh = fit$llh
    12. 

  12.     .params$params[INDEX] = fit$par
    13. 

  13.     .setfGarchEnv(.params = .params)
    14. 

  14.   }
    15. 

  15.   .params$llh = fit$llh
    16. 

  16.   .params$params[INDEX] = fit$par
    17. 

  17.   .setfGarchEnv(.params = .params)
    18. 

  18.   if (hessian == "ropt") {
    19. 

  19.     fit$hessian <- -.garchRoptimhess(par = fit$par, .params = .params, 
    20. 

  20.       .series = .series)
    21. 

  21.     titleHessian = "R-optimhess"
    22. 

  22.   }
    23. 

  23. ...............
    24. 

  24.   if (.params$control$xscale) {
    25. 

  25.     .series$x <- .series$x * .series$scale          #将标准化后的值还原为原始的数值了
    26. 

  26.     if (.params$include["mu"]) 
    27. 

  27.       fit$coef["mu"] <- fit$par["mu"] <- .params$params["mu"] <- .params$params["mu"] * 
    28. 

  28.         .series$scale
    29. 

  29.     if (.params$include["omega"]) 
    30. 

  30.       fit$coef["omega"] <- fit$par["omega"] <- .params$params["omega"] <- .params$params["omega"] * 
    31. 

  31.         .series$scale^(.params$params["delta"])
    32. 

  32.     .setfGarchEnv(.params = .params)
    33. 

  33.     .setfGarchEnv(.series = .series)
    34. 

  34.   }
    35. 

  35.   if (.params$control$xscale) {
    36. 

  36.     if (.params$include["mu"]) {
    37. 

  37.       fit$hessian[, "mu"] <- fit$hessian[, "mu"]/.series$scale
    38. 

  38.       fit$hessian["mu", ] <- fit$hessian["mu", ]/.series$scale
    39. 

  39.     }
    40. 

  40.     if (.params$include["omega"]) {
    41. 

  41.       fit$hessian[, "omega"] <- fit$hessian[, "omega"]/.series$scale^(.params$params["delta"])
    42. 

  42.       fit$hessian["omega", ] <- fit$hessian["omega", ]/.series$scale^(.params$params["delta"])
    43. 

  43.     }
    44. 

  44.   }
    45. 

  45.   .llh <- fit$llh <- fit$value <- .garchLLH(fit$par, trace = FALSE, 
    46. 

  46.     fGarchEnv = TRUE)
    47. 

  47.   .series <- .getfGarchEnv(".series")
    48. 

  48.   if (robust.cvar) 
    49. 

  49.     fit$gradient <- -.garchRCDAGradient(par = fit$par, .params = .params, 
    50. 

  50.       .series = .series)
    51. 

  51.   N = length(.series$x)
    52. 

  52.   NPAR = length(fit$par)
    53. 

  53.   fit$ics = c(AIC = c((2 * fit$value)/N + 2 * NPAR/N), BIC = (2 * 
    54. 

  54.     fit$value)/N + NPAR * log(N)/N, SIC = (2 * fit$value)/N + 
    55. 

  55.     log((N + 2 * NPAR)/N), HQIC = (2 * fit$value)/N + (2 * 
    56. 

  56.     NPAR * log(log(N)))/N)
    57. 

  57.   names(fit$ics) <- c("AIC", "BIC", "SIC", "HQIC")
    58. 

  58.   fit
    59. 

  59. }
    60.

其中4 (4.0.2.4.1).garchRnlminb函数
.garchRnlminb函数
    

  1. function (.params, .series, .garchLLH, trace) 
    2. 

  2. {
    3. 

  3.   if (trace) 
    4. 

  4.     cat("\nR coded nlminb Solver: \n\n")
    5. 

  5.   INDEX = .params$index
    6. 

  6.   parscale = rep(1, length = length(INDEX))
    7. 

  7.   names(parscale) = names(.params$params[INDEX])
    8. 

  8.   parscale["omega"] = var(.series$x)^(.params$delta/2)
    9. 

  9.   parscale["mu"] = abs(mean(.series$x))           #这个并非mu的初始值
    10. 

  10.   TOL1 = .params$control$tol1
    11. 

  11.   fit <- nlminb(start = .params$params[INDEX], objective = .garchLLH, 
    12. 

  12.     lower = .params$U[INDEX], upper = .params$V[INDEX], 
    13. 

  13.     scale = 1/parscale, control = list(eval.max = 2000, 
    14. 

  14.       iter.max = 1500, rel.tol = 1e-14 * TOL1, x.tol = 1e-14 * 
    15. 

  15.         TOL1, trace = as.integer(trace)), fGarchEnv = FALSE)
    16. 

  16.   fit$value <- fit.llh <- fit$objective
    17. 

  17.   names(fit$par) = names(.params$params[INDEX])
    18. 

  18.   .garchLLH(fit$par, trace = FALSE, fGarchEnv = TRUE)
    19. 

  19.   fit$coef <- fit$par
    20. 

  20.   fit$llh <- fit$objective
    21. 

  21.   fit
    22. 

  22. }
    23.
关于nlminb函数，请看笔记：
.params$params[INDEX]是初始值
 objective = .garchLLH是被优化的函数
数据呢？怎么没见传入数据？

.garchRnlminb函数函数运行完，则会估计出参数：
fit$par
        mu      omega     alpha1      beta1 
0.08876900 0.05705989 0.08643831 0.85258554 
这个值呢，是别标准化后的参数值
  1. 
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函数
    

  1. .garchLLH
    2. 

  2. function (params, trace = TRUE, fGarchEnv = FALSE) 
    3. 

  3. {
    4. 

  4.     DEBUG = FALSE
    5. 

  5.     if (DEBUG) 
    6. 

  6.         print("Entering Function .garchLLH")
    7. 

  7.     .series <- .getfGarchEnv(".series")
    8. 

  8.     .params <- .getfGarchEnv(".params")
    9. 

  9.     .garchDist <- .getfGarchEnv(".garchDist")   #概率密度函数
    10. 

  10.     .llh <- .getfGarchEnv(".llh")  #对数似然函数
    11. 

  11.     if (DEBUG) 
    12. 

  12.         print(.params$control$llh)
    13. 

  13.     if (.params$control$llh == "internal") {
    14. 

  14.         if (DEBUG) 
    15. 

  15.             print("internal")
    16. 

  16.         return(.aparchLLH.internal(params, trace = trace, fGarchEnv = fGarchEnv))  #结束
    17. 

  17.     }
    18. 

  18.     else if (.params$control$llh == "filter") {
    19. 

  19.         if (DEBUG) 
    20. 

  20.             print("filter")
    21. 

  21.         return(.aparchLLH.filter(params, trace = trace, fGarchEnv = fGarchEnv))
    22. 

  22.     }
    23. 

  23.     else if (.params$control$llh == "testing") {
    24. 

  24.         if (DEBUG) 
    25. 

  25.             print("testing")
    26. 

  26.         return(.aparchLLH.testing(params, trace = trace, fGarchEnv = fGarchEnv))
    27. 

  27.     }
    28. 

  28.     else {
    29. 

  29.         stop("LLH is neither internal, testing, nor filter!")
    30. 

  30.     }
    31. 

  31. }
    32. 

  32. <environment: namespace:fGarch>
    33.
.aparchLLH.internal函数
其中我进去看了下
    

  1. function (params, trace = TRUE, fGarchEnv = TRUE) 
    2. 

  2. {
    3. 

  3.     DEBUG = FALSE
    4. 

  4.     if (DEBUG) 
    5. 

  5.         print("Entering Function .garchLLH.internal")
    6. 

  6.     .series <- .getfGarchEnv(".series")
    7. 

  7.     .params <- .getfGarchEnv(".params")
    8. 

  8.     .garchDist <- .getfGarchEnv(".garchDist")
    9. 

  9.     .llh <- .getfGarchEnv(".llh")
    10. 

  10.     if (DEBUG) 
    11. 

  11.         print(.params$control$llh)
    12. 

  12.     if (.params$control$llh == "internal") {
    13. 

  13.         INDEX <- .params$index
    14. 

  14.         MDIST <- c(norm = 10, QMLE = 10, snorm = 11, std = 20, 
    15. 

  15.             sstd = 21, ged = 30, sged = 31)[.params$cond.dist]
    16. 

  16.         if (.params$control$fscale) 
    17. 

  17.             NORM <- length(.series$x)
    18. 

  18.         else NORM = 1
    19. 

  19.         REC <- 1
    20. 

  20.         if (.series$init.rec == "uev") 
    21. 

  21.             REC <- 2
    22. 

  22.         MYPAR <- c(REC = REC, LEV = as.integer(.params$leverage), 
    23. 

  23.             MEAN = as.integer(.params$includes["mu"]), DELTA = as.integer(.params$includes["delta"]), 
    24. 

  24.             SKEW = as.integer(.params$includes["skew"]), SHAPE = as.integer(.params$includes["shape"]), 
    25. 

  25.             ORDER = .series$order, NORM = as.integer(NORM))
    26. 

  26.         MAX <- max(.series$order)
    27. 

  27.         NF <- length(INDEX)
    28. 

  28.         N <- length(.series$x)
    29. 

  29.         DPARM <- c(.params$delta, .params$skew, .params$shape)
    30. 

  30.         fit <- .Fortran("garchllh", N = as.integer(N), Y = as.double(.series$x), 
    31. 

  31.             Z = as.double(.series$z), H = as.double(.series$h), 
    32. 

  32.             NF = as.integer(NF), X = as.double(params), DPARM = as.double(DPARM), 
    33. 

  33.             MDIST = as.integer(MDIST), MYPAR = as.integer(MYPAR), 
    34. 

  34.             F = as.double(0), PACKAGE = "fGarch")
    35. 

  35.         llh <- fit[[10]]                            #此时llh就是-312
    36. 

  36.         if (is.na(llh)) 
    37. 

  37.             llh = .llh + 0.1 * (abs(.llh))
    38. 

  38.         if (!is.finite(llh)) 
    39. 

  39.             llh = .llh + 0.1 * (abs(.llh))
    40. 

  40.         .setfGarchEnv(.llh = llh)
    41. 

  41.         if (fGarchEnv) {
    42. 

  42.             .series$h <- fit[[4]]
    43. 

  43.             .series$z <- fit[[3]]                   #直到这里，z才从444个0被替换为残差
    44. 

  44.             .setfGarchEnv(.series = .series)
    45. 

  45.         }
    46. 

  46.     }
    47. 

  47.     else {
    48. 

  48.         stop("LLH is not internal!")
    49. 

  49.     }
    50. 

  50.     c(LogLikelihood = llh)
    51. 

  51. }
    52.
.Fortran是调用Fortran代码的，真是我擦！

    

  1. Browse[5]> .garchDist
    2. 

  2. function (z, hh, skew, shape) 
    3. 

  3. {
    4. 

  4.     dnorm(x = z/hh, mean = 0, sd = 1)/hh
    5. 

  5. }
    6. 

  6. <environment: namespace:fGarch>
    7.


写在2017年1月28日00:54:55
农历新年凌晨，新的一年要更加努力~

使用RStudio调试(debug)基础学习(二)和fGarch包中的garchFit函数估计GARCH模型的原理和源码的更多相关文章

  1. 使用RStudio调试&lpar;debug&rpar;基础学习&lpar;一&rpar;

    点击行号的左侧,即可设置断点(或者按下Shift+F9),如果没有出现,反而出现下图的警告: 那么只是因为我的坏习惯--写一段脚本测试的时候都是新建,但不save到本地,不喜欢保存,写的差不多了才开始 ...

  2. Python入门基础学习 二

    Python入门基础学习 二 猜数字小游戏进阶版 修改建议: 猜错的时候程序可以给出提示,告诉用户猜测的数字偏大还是偏小: 没运行一次程序只能猜测一次,应该提供多次机会给用户猜测: 每次运行程序,答案 ...

  3. Python基础学习二

    Python基础学习二 1.编码 utf-8编码:自动将英文保存为1个字符,中文3个字符.ASCll编码被囊括在内. unicode:将所有字符保存为2给字符,容纳了世界上所有的编码. 2.字符串内置 ...

  4. Go基础学习&lpar;二&rpar;

    数组[array] 数组定义[定义后长度不可变] 12 symbol := [...]string{USD: "$", EUR: "€", GBP: &quot ...

  5. Django基础学习二

    今天继续学习django的基础 学习用户提交url如何获得返回值 1.首先需要在工程的urls文件定义指定的urls要路由给哪个函数 在这个例子中,我们定义home的urls路由给views里的tes ...

  6. salesforce lightning零基础学习&lpar;二&rpar; lightning 知识简单介绍----lightning事件驱动模型

    看此篇博客前或者后,看一下trailhead可以加深印象以及理解的更好:https://trailhead.salesforce.com/modules/lex_dev_lc_basics 做过cla ...

  7. CSS入门基础学习二

    我们下午继续学习CSS的入门基础,搬上你的小板凳赶快进入吧! 一.背景(background) Background-color:背景颜色 background-image (背景图片) backgr ...

  8. WebService基础学习&lpar;二&rpar;&mdash&semi;三要素

    一.Java中WebService规范      JAVA *有三种WebService 规范,分别是JAX-WS.JAX-RS.JAXM&SAAJ(废弃).   1.JAX-WS规范    ...

  9. jQuery基础学习&lpar;二&rpar;&mdash&semi;jQuery选择器

    一.jQuery基本选择器 1.CSS选择器     在学习jQuery选择器之前,先介绍一下之前学过的CSS选择器. 选择器 语法 描述 示例   标签选择器 E {                 ...

随机推荐

  1. 6&period;DNS公司PC访问外网的设置 &plus; 主DNS服务器和辅助DNS服务器的配置

    网站部署之~Windows Server | 本地部署 http://www.cnblogs.com/dunitian/p/4822808.html#iis DNS服务器部署不清楚的可以看上一篇:ht ...

  2. JS函数节流

    背景:在前端开发中,有时会为页面绑定resize事件,或为一个页面元素拖拽事件(其核心就是绑定mousemove)在一个正常操作中也有可能在一个短时间内触发非常多次事件绑定程序,而DOM操作是很消耗性 ...

  3. 有关uploadifive的使用经验

    这段时间做了一个项目用到uploadifive上传控件,和uploadify不同,这个控件是基于HTML5的版本而不用支持falsh,因而移动端也可以使用. 整理用过的相关属性与方法: 属性 作用 a ...

  4. XE3随想14：关于 SO 与 SA 函数

    通过 SuperObject 的公用函数 SO 实现一个 ISuperObject 接口非常方便; 前面都是给它一个字符串参数, 它的参数可以是任一类型甚至是常数数组. SA 和 SO 都是返回一 I ...

  5. 损失函数&lpar;Loss Function&rpar; -1

    http://www.ics.uci.edu/~dramanan/teaching/ics273a_winter08/lectures/lecture14.pdf Loss Function 损失函数 ...

  6. Flask-WTF form doesn&&num;39&semi;t have attribute &&num;39&semi;validate&lowbar;on&lowbar;submit&&num;39&semi;问题

    今天在学习WTF表单的时候遇到了这个问题,在*上搜索查到了解决方案 from flask.ext.wtf import Form from wtforms import Tex ...

  7. 文件上传时jquery&period;form&period;js中提示form&period;submit SCRIPT5&colon; 拒绝访问

    利用其它控件触发file的click事件来选择文件后,使用jquery.form.js中的submit方法提交时IE报错:form.submit SCRIPT5: 拒绝访问,其它浏览器正常, < ...

  8. Session（）

    如何使用 Session Java Api 只给我们一种方式来 获取 当前会话相关的 session: HttpSession session = request.getSession(); //或 ...

  9. Web Worker 初探

    什么是Web Worker? Web Worker 是Html5 提出的能够在后台运行javascript的对象,独立于其他脚本,不会影响页面的性能,也不会影响你继续对于页面进行操作.通俗点讲,就是后 ...

  10. 移动端js触摸事件大全

    一.手机上的触摸事件 基本事件: touchstart   //手指刚接触屏幕时触发 touchmove    //手指在屏幕上移动时触发 touchend     //手指从屏幕上移开时触发 下面这 ...

相关文章