（持续更新）学习R语言那些事

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

此贴用于记录学习R语言记录

1.普通多元回归分析的R语言代码

#设置工作路径
setwd("C:/Users/smile/Desktop/yoga/20220419")

#安装必须包
install.packages("openxlsx")
install.packages("ggcorrplot")
install.packages("ggthemes")
install.packages("car")
install.packages("stargazer")

#读取安装包
library(openxlsx)
library(ggcorrplot)
library(ggthemes)
library(car)
library(stargazer)

#读取数据集
data <- read.xlsx("发送给老师.xlsx",sheet = 1)
data <- na.omit(data)

#1.描述性统计
summary(data)

#2.相关性分析
#矩阵+热图
corr <- round(cor(data[,c(3:11)],method = c("spearman")),3)

pdf("相关性图像.pdf")
ggcorrplot(corr,method = "square",hc.method = "ward.D",
           outline.color = "white",ggtheme = theme_bw(),
           type = "upper",colors = c("blue","white","red"),
           lab = T,lab_size = 5)
dev.off()

#导出相关性矩阵
write.csv(corr,"相关性系数矩阵.csv")

#3.逐步回归
model <- lm(y~.,data = data[,c(3:11)])
step_model <- step(model)
summary(step_model)

#4.共线性分析
vif(step_model)
#最终的特征VIF均小于10，说明x1,x3,x5,x6与y均不存在严重共线性。

#5.补齐所有逐步回归的过程数据
model1 <- lm(y~ xx + x1 + x2 + x3 + x4 + x5 + x6,data = data[,c(3:11)])
model2 <- lm(y~ xx + x1 + x2 + x3 + x5 + x6,data = data[,c(3:11)])
model3 <- lm(y~ x1 + x2 + x3 + x5 + x6,data = data[,c(3:11)])

#导出回归模型
stargazer( model,model1,model2,model3,step_model,
           title = "results", align = F, type = "text",
           no.space = TRUE, out = "fit.html")

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：R语言 那些事 R语言学习

2 固定效应回归的R语言代码

3.学习R遇到一些问题
could not find function "kable"

2.R 语言  GARCH族代码分析
代码：
getwd()
mydata=read.table('r.txt',header=T)
library(View)
View(mydata)
mydata$time=as.Date(mydata$time)


install.packages('PerformanceAnalytics')
install.packages('CalculateReturns')
install.packages('forecast')


library(xts)
library(PerformanceAnalytics)
library(CalculateReturns)
library(urca)
library(fGarch)
library(forecast)
library(copula)
library(tseries)

mydata=xts(mydata[,-1],order.by = mydata$time)

rt=coredata(mydata['1986-01-03/2021-12-31']$rt)

date=index(mydata)
l=dim(mydata)[1]
### plot returns
par(mfcol=c(2,2),mar=c(2,2,2,2))
r=coredata(mydata['1986-01-03/2021-12-31']$r)

g1<-garch(r, order = c(1, 1))
summary(g1)


md_name <- c("sGARCH", "eGARCH", "gjrGARCH",
             "apARCH", "TGARCH", )
                         
md_dist <- c("norm", "snorm","std", "sstd")


# 设定函数
model_train <- function(data_, m_name, m_dist) {
  res <- NULL
  for (i in m_name) {
    sub_model <- NULL
    main_model <- i
    if (i == "TGARCH" | i == "NAGARCH") {
      sub_model <- i
      main_model <- "fGARCH"
    }

    for (j in m_dist) {
      model_name <- paste0('ARMA(4,4)', '-',
                           i, '-', j,
                           sep = "")
      if (main_model == "fGARCH") {
        model_name <- paste0('ARMA(4,4)', '-',
                             sub_model, '-', j,
                             sep = "")
      }
      print(model_name)
      # 模型设定
      mean.spec <- list(armaOrder = c(4, 4), include.mean = F,
                        archm = F, archpow = 1, arfima = F,
                        external.regressors = NULL)
      var.spec <- list(model = main_model, garchOrder = c(1, 1),
                       submodel = sub_model,
                       external.regressors = NULL,
                       variance.targeting = F)
      dist.spec <- j
      my_spec <- ugarchspec(mean.model = mean.spec,
                            variance.model = var.spec,
                            distribution.model = dist.spec)
      # 模型拟合
      my_fit <- ugarchfit(data = data_, spec = my_spec)
      res <- rbind(res, c(list(ModelName = model_name),
                          list(LogL = likelihood(my_fit)),
                          list(AIC = infocriteria(my_fit)[1]),
                          list(BIC = infocriteria(my_fit)[2]),
                          list(RMSE = sqrt(mean(residuals(my_fit)^2))),
                          list(md = my_fit)))
    }
  }
  return(res)
}

# 模型拟合
my_model <- model_train(r.data, md_name, md_dist)

# 保存模型结果
write.table(my_model[, 1:5],
            "D:\\Rproject\\result.txt",
            sep = " ")

# 读取模型结果
res_table <- read.table("D:\\Rproject\\result.txt",
                        header = TRUE,
                        sep = " ")
kable(res_table, caption = "Result of GARCH Model", align = 'l')

# 选择模型再拟合
md_name <- c("eGARCH")
md_dist <- c("norm", "snorm", "std", "sstd", "ged", "sged", "jsu")
new_model <- model_train(r.data, md_name, md_dist)

# 信息准则表
my_info_cri <- NULL
for (i in 1:dim(new_model)[1]) {
  new_aic <- round(infocriteria(new_model[i, 6]$md)[1], digits = 4)
  new_bic <- round(infocriteria(new_model[i, 6]$md)[2], digits = 4)
  new_sib <- round(infocriteria(new_model[i, 6]$md)[3], digits = 4)
  new_hq <- round(infocriteria(new_model[i, 6]$md)[4], digits = 4)
  new_llk <- round(likelihood(new_model[i, 6]$md), digits = 4)
  my_info_cri <- rbind(my_info_cri, c(new_model[i, 1],
                                      Akaike = new_aic,
                                      Bayes = new_bic,
                                      Shibata = new_sib,
                                      HQ = new_hq,
                                      LLH = new_llk))
}


# 模型系数显著性表
my_cof_table <- NULL
my_cof_names <- NULL
for (i in 1:dim(new_model)[1]) {
  my_cof_names <- rbind(my_cof_names, new_model[i, 1]$ModelName)
  my_df <- data.frame(t(round(new_model[i, 6]$md@fit$robust.matcoef[, 4],
                              digits = 4)))
  if (i == 1) {
    my_cof_table <- data.frame(my_df)
  }
  else {
    my_cof_table <- full_join(my_cof_table, my_df)
  }
}
row.names(my_cof_table) <- my_cof_names[, 1]
my_cof_table <- t(my_cof_table)


# 参数稳定性个别检验表
my_nyblom_table <- NULL
my_nyblom_name <- NULL
for (i in 1:dim(new_model)[1]) {
  my_nyblom_name <- rbind(my_nyblom_name,
                          new_model[i, 1]$ModelName)
  my_df <- data.frame(t(round(nyblom(new_model[i, 6]$md)$IndividualStat,
                              digits = 4)))
  if (i == 1) {
    my_nyblom_table <- data.frame(my_df)
  }
  else {
    my_nyblom_table <- full_join(my_nyblom_table, my_df)
  }
}
row.names(my_nyblom_table) <- my_nyblom_name[, 1]
my_IC <- NULL
for (i in 1:dim(new_model)[1]) {
  my_IC <- cbind(my_IC, nyblom(new_model[1, 6]$md)$IndividualCritical)
}
my_nyblom_table <- rbind(t(my_nyblom_table), my_IC)

# 参数稳定性联合检验表
my_nyblom_table2 <- NULL
my_nyblom_name2 <- NULL
for (i in 1:dim(new_model)[1]) {
  my_nyblom_name2 <- rbind(my_nyblom_name2,
                           new_model[i, 1]$ModelName)
  df1 <- data.frame(JoinStat = round(nyblom(new_model[i,6]$md)$JointStat,
                                     digits = 4))
  df2 <- data.frame(mm = round(nyblom(new_model[i, 6]$md)$JointCritical,
                               digits = 4))
  temp_table <- cbind(df1, t(df2))
  if (i == 1) {
    my_nyblom_table2 <- temp_table
  }
  else {
    my_nyblom_table2 <- full_join(my_nyblom_table2, temp_table)
  }
}
row.names(my_nyblom_table2) <- my_nyblom_name2

# 模型分布拟合度p值检验表
my_gof_table <- NULL
my_gof_name <- NULL
group_name <- c(20, 30, 40, 50)
for (i in 1:dim(new_model)[1]) {
  my_gof_name <- rbind(my_gof_name,
                       new_model[i, 1]$ModelName)
  df1 <- round(gof(new_model[i, 6]$md,
                   groups = group_name)[1:4, 3],
               digits = 4)
  if (i == 1) {
    my_gof_table <- df1
  }
  else {
    my_gof_table <- rbind(my_gof_table, df1)
  }
}
row.names(my_gof_table) <- my_gof_name
colnames(my_gof_table) <- c(20, 30, 40, 50)

# 符号偏误显著性检验表
my_sign_table <- NULL
my_sign_name <- NULL
row_sign_name <- row.names(signbias(new_model[1, 6]$md))
for (i in 1:dim(new_model)[1]) {
  my_sign_name <- rbind(my_sign_name,
                        new_model[i, 1]$ModelName)
  df1 <- round(signbias(new_model[i, 6]$md)$prob, digits = 4)
  if (i == 1) {
    my_sign_table <- df1
  }
  else {
    my_sign_table <- rbind(my_sign_table, df1)
  }
}
row.names(my_sign_table) <- my_sign_name
colnames(my_sign_table) <- row_sign_name

# 展示表格
kable(my_cof_table,
      caption = "P-value Table of Coefficients")
kable(my_nyblom_table,
      caption = "P-value Table of Nyblom Individual Stability Test")
kable(my_nyblom_table2,
      caption = "P-value Table of Nyblom Joint Stability Test")
kable(my_sign_table,
      caption = "P-value Table of Sign Bias Test")
kable(my_gof_table,
      caption = "P-value Table of Goodness-of-fit")
kable(my_info_cri,
      caption = "Table of Information Criteria")

# 最终模型系数表
kable(new_model[6, 6]$md@fit$robust.matcoef,
      caption = "Table of Final Model Coef.")

所有的代码都是亲测跑过的，需要答疑  可以私我

md_name <- c("sGARCH", "eGARCH", "gjrGARCH",
             "apARCH", "TGARCH", )
                         
md_dist <- c("norm", "snorm","std", "sstd")


# 设定函数
model_train <- function(data_, m_name, m_dist) {
  res <- NULL
  for (i in m_name) {
    sub_model <- NULL
    main_model <- i
    if (i == "TGARCH" | i == "NAGARCH") {
      sub_model <- i
      main_model <- "fGARCH"
    }

    for (j in m_dist) {
      model_name <- paste0('ARMA(4,4)', '-',
                           i, '-', j,
                           sep = "")
      if (main_model == "fGARCH") {
        model_name <- paste0('ARMA(4,4)', '-',
                             sub_model, '-', j,
                             sep = "")
      }
      print(model_name)
      # 模型设定
      mean.spec <- list(armaOrder = c(4, 4), include.mean = F,
                        archm = F, archpow = 1, arfima = F,
                        external.regressors = NULL)
      var.spec <- list(model = main_model, garchOrder = c(1, 1),
                       submodel = sub_model,
                       external.regressors = NULL,
                       variance.targeting = F)
      dist.spec <- j
      my_spec <- ugarchspec(mean.model = mean.spec,
                            variance.model = var.spec,
                            distribution.model = dist.spec)
      # 模型拟合
      my_fit <- ugarchfit(data = data_, spec = my_spec)
      res <- rbind(res, c(list(ModelName = model_name),
                          list(LogL = likelihood(my_fit)),
                          list(AIC = infocriteria(my_fit)[1]),
                          list(BIC = infocriteria(my_fit)[2]),
                          list(RMSE = sqrt(mean(residuals(my_fit)^2))),
                          list(md = my_fit)))
    }
  }
  return(res)
}

# 模型拟合
my_model <- model_train(r.data, md_name, md_dist)

# 保存模型结果
write.table(my_model[, 1:5],
            "D:\\Rproject\\result.txt",
            sep = " ")

# 读取模型结果
res_table <- read.table("D:\\Rproject\\result.txt",
                        header = TRUE,
                        sep = " ")
kable(res_table, caption = "Result of GARCH Model", align = 'l')

# 选择模型再拟合
md_name <- c("eGARCH")
md_dist <- c("norm", "snorm", "std", "sstd", "ged", "sged", "jsu")
new_model <- model_train(r.data, md_name, md_dist)

# 信息准则表
my_info_cri <- NULL
for (i in 1:dim(new_model)[1]) {
  new_aic <- round(infocriteria(new_model[i, 6]$md)[1], digits = 4)
  new_bic <- round(infocriteria(new_model[i, 6]$md)[2], digits = 4)
  new_sib <- round(infocriteria(new_model[i, 6]$md)[3], digits = 4)
  new_hq <- round(infocriteria(new_model[i, 6]$md)[4], digits = 4)
  new_llk <- round(likelihood(new_model[i, 6]$md), digits = 4)
  my_info_cri <- rbind(my_info_cri, c(new_model[i, 1],
                                      Akaike = new_aic,
                                      Bayes = new_bic,
                                      Shibata = new_sib,
                                      HQ = new_hq,
                                      LLH = new_llk))
}


# 模型系数显著性表
my_cof_table <- NULL
my_cof_names <- NULL
for (i in 1:dim(new_model)[1]) {
  my_cof_names <- rbind(my_cof_names, new_model[i, 1]$ModelName)
  my_df <- data.frame(t(round(new_model[i, 6]$md@fit$robust.matcoef[, 4],
                              digits = 4)))
  if (i == 1) {
    my_cof_table <- data.frame(my_df)
  }
  else {
    my_cof_table <- full_join(my_cof_table, my_df)
  }
}
row.names(my_cof_table) <- my_cof_names[, 1]
my_cof_table <- t(my_cof_table)


# 参数稳定性个别检验表
my_nyblom_table <- NULL
my_nyblom_name <- NULL
for (i in 1:dim(new_model)[1]) {
  my_nyblom_name <- rbind(my_nyblom_name,
                          new_model[i, 1]$ModelName)
  my_df <- data.frame(t(round(nyblom(new_model[i, 6]$md)$IndividualStat,
                              digits = 4)))
  if (i == 1) {
    my_nyblom_table <- data.frame(my_df)
  }
  else {
    my_nyblom_table <- full_join(my_nyblom_table, my_df)
  }
}
row.names(my_nyblom_table) <- my_nyblom_name[, 1]
my_IC <- NULL
for (i in 1:dim(new_model)[1]) {
  my_IC <- cbind(my_IC, nyblom(new_model[1, 6]$md)$IndividualCritical)
}
my_nyblom_table <- rbind(t(my_nyblom_table), my_IC)

# 参数稳定性联合检验表
my_nyblom_table2 <- NULL
my_nyblom_name2 <- NULL
for (i in 1:dim(new_model)[1]) {
  my_nyblom_name2 <- rbind(my_nyblom_name2,
                           new_model[i, 1]$ModelName)
  df1 <- data.frame(JoinStat = round(nyblom(new_model[i,6]$md)$JointStat,
                                     digits = 4))
  df2 <- data.frame(mm = round(nyblom(new_model[i, 6]$md)$JointCritical,
                               digits = 4))
  temp_table <- cbind(df1, t(df2))
  if (i == 1) {
    my_nyblom_table2 <- temp_table
  }
  else {
    my_nyblom_table2 <- full_join(my_nyblom_table2, temp_table)
  }
}
row.names(my_nyblom_table2) <- my_nyblom_name2

# 模型分布拟合度p值检验表
my_gof_table <- NULL
my_gof_name <- NULL
group_name <- c(20, 30, 40, 50)
for (i in 1:dim(new_model)[1]) {
  my_gof_name <- rbind(my_gof_name,
                       new_model[i, 1]$ModelName)
  df1 <- round(gof(new_model[i, 6]$md,
                   groups = group_name)[1:4, 3],
               digits = 4)
  if (i == 1) {
    my_gof_table <- df1
  }
  else {
    my_gof_table <- rbind(my_gof_table, df1)
  }
}
row.names(my_gof_table) <- my_gof_name
colnames(my_gof_table) <- c(20, 30, 40, 50)

# 符号偏误显著性检验表
my_sign_table <- NULL
my_sign_name <- NULL
row_sign_name <- row.names(signbias(new_model[1, 6]$md))
for (i in 1:dim(new_model)[1]) {
  my_sign_name <- rbind(my_sign_name,
                        new_model[i, 1]$ModelName)
  df1 <- round(signbias(new_model[i, 6]$md)$prob, digits = 4)
  if (i == 1) {
    my_sign_table <- df1
  }
  else {
    my_sign_table <- rbind(my_sign_table, df1)
  }
}
row.names(my_sign_table) <- my_sign_name
colnames(my_sign_table) <- row_sign_name

# 展示表格
kable(my_cof_table,
      caption = "P-value Table of Coefficients")
kable(my_nyblom_table,
      caption = "P-value Table of Nyblom Individual Stability Test")
kable(my_nyblom_table2,
      caption = "P-value Table of Nyblom Joint Stability Test")
kable(my_sign_table,
      caption = "P-value Table of Sign Bias Test")
kable(my_gof_table,
      caption = "P-value Table of Goodness-of-fit")
kable(my_info_cri,
      caption = "Table of Information Criteria")

# 最终模型系数表
kable(new_model[6, 6]$md@fit$robust.matcoef,
      caption = "Table of Final Model Coef.")

久菜盒子工作室 发表于 2022-8-6 16:33
3.学习R遇到一些问题
could not find function "kable"

请问这个问题怎么解决？