华裔计量经济学家 Changyou Sun 教授的 erer 工具包推荐

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

欢迎使用Markdown编辑器

经管之家：Do the best economic and management education！
Author @ bbs.pinggu.org: Daniel Tulpen Liu ！
Author @ RStudio | Open source & professional software for data science teams - RStudio : Changyou Sun (2020). erer: Empirical Research in Economics with R.
Copyright © Changyou Sun (2020)

你好！ 这是你第一次使用 Markdown编辑器 所展示的欢迎页。如果你想学习如何使用Markdown编辑器, 可以仔细阅读这篇文章，了解一下Markdown的基本语法知识。

今天要分享的是 MATLAB 的程序；

作者是国际上特别牛的人， 程序的实现的话， 我认为没有十年的 MATLAB 使用经验的话，是做不到的。
所以呢，只有好好的学习。这个代码的分享是为了大家的学习。

• 复杂的是，这个MATLAB 程序还调用了 JAVA 程序，可见作者的编程功底非常的霸道。
• Enjoy

Mani routine of the package

#  Author :
#  Title  :
#  E-mail :


# References 'packages' ---------------------------------------------------

citation("fBasics")
citation("urca")

citation("erer")

# Initial Envirenment -----------------------------------------------------
rm(list = ls(all=TRUE))
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

opar = par(mfrow = c(1, 1))


# load package ------------------------------------------------------------

library(erer)

# one ---------------------------------------------------------------------
# A new function: probabilities for ordered choice
ocProb <- function(w, nam.c, n = 100, digits = 3)
{
  # 1. Check inputs
  if (!inherits(w, "polr")) {
    stop("Need an ordered choice model from 'polr()'.\n")
  }
  if (w$method != "probit" & w$method != "logistic") {
    stop("Need a probit or logit model.\n")
  }
  if (missing(nam.c)) stop("Need a continous variable name'.\n")
  
  # 2. Abstract data out
  lev <- w$lev; J <- length(lev)
  x.name <- attr(x = w$terms, which = "term.labels")
  x2 <- w$model[, x.name]
  if (identical(sort(unique(x2[, nam.c])), c(0, 1)) ||
      inherits(x2[, nam.c], what = "factor")) {
    stop("nam.c must be a continuous variable.")
  }
  
  ww <- paste("~ 1", paste("+", x.name, collapse = " "), collapse = " ")
  x <- model.matrix(as.formula(ww), data = x2)[, -1]
  b.est <- as.matrix(coef(w)); K <- nrow(b.est)
  z <- c(-10^6, w$zeta, 10^6)  # expand it with two extreme thresholds
  z2 <- matrix(data = z, nrow = n, ncol = length(z), byrow = TRUE)
  
  pfun <- switch(w$method, probit = pnorm, logistic = plogis)
  dfun <- switch(w$method, probit = dnorm, logistic = dlogis)
  V2 <- vcov(w)  # increase covarance matrix by 2 fixed thresholds
  V3 <- rbind(cbind(V2, 0, 0), 0, 0)
  ind <- c(1:K, nrow(V3)-1, (K+1):(K+J-1), nrow(V3))
  V4 <- V3[ind, ]; V5 <- V4[, ind]
  
  # 3. Construct x matrix and compute xb
  mm <- matrix(data = colMeans(x), ncol = ncol(x), nrow = n, byrow = TRUE)
  colnames(mm) <- colnames(x)
  ran <- range(x[, nam.c])
  mm[, nam.c] <- seq(from = ran[1], to = ran[2], length.out = n)
  xb <- mm %*% b.est
  xb2 <- matrix(data = xb, nrow = n, ncol = J, byrow = FALSE)  # J copy
  
  # 4. Compute probability by category; vectorized on z2 and xb2
  pp <- pfun(z2[, 2:(J+1)] - xb2) - pfun(z2[, 1:J] - xb2)
  trend <- cbind(mm[, nam.c], pp)
  colnames(trend) <- c(nam.c, paste("p", lev, sep="."))
  
  # 5. Compute the standard errors
  se <- matrix(data = 0, nrow = n, ncol = J)
  for (i in 1:J) {
    z1 <- z[i] - xb; z2 <- z[i+1] - xb
    d1 <- diag(c(dfun(z1) - dfun(z2)), n, n) %*% mm
    q1 <- - dfun(z1); q2 <-   dfun(z2)
    dr <- cbind(d1, q1, q2)
    V <- V5[c(1:K, K+i, K+i+1), c(1:K, K+i, K+i+1)]
    va <- dr %*% V %*% t(dr)
    se[, i] <- sqrt(diag(va))
  }
  colnames(se) <- paste("Pred_SE", lev, sep = ".")
  
  # 6. Report results
  t.value <- pp / se
  p.value <- 2 * (1 - pt(abs(t.value), n - K))
  out <- list()
  for (i in 1:J) {
    out[[i]] <- round(cbind(predicted_prob = pp[, i], error = se[, i],
                            t.value = t.value[, i], p.value = p.value[, i]), digits)
  }
  out[[J+1]] <- round(x = trend, digits = digits)
  names(out) <- paste("predicted_prob", c(lev, "all"), sep = ".")
  result <- listn(w, nam.c, method=w$method, mean.x=colMeans(x), out, lev)
  class(result) <- "ocProb"; return(result)
}

# Example: include "Freq" to have a continuous variable for demo
library(erer); library(MASS); data(housing); str(housing); tail(housing)
reg2 <- polr(formula = Sat ~ Infl + Type + Cont + Freq, data = housing, 
             Hess = TRUE, method = "probit")
p2 <- ocProb(w = reg2, nam.c = 'Freq', n = 300); p2
plot(p2)


# two ---------------------------------------------------------------------

# Shared commands for versions A and B
curve(expr = -k - b * sqrt(1 - ((x - h) / a)^2), 
      from = 0, to = 10, n = 500, add = TRUE) 
box()
segments(x0 = 0, y0 = 0, x1 = 10, y1 = 0, col = "gray80", lwd = 10, lend=2)
arrows(x0 = 1, y0 = c(-1.5, 1.5), x1 = 1, y1 = c(-3,3), col = "gray80",
       lwd = 10, code = 2, length = 0.15, angle = 20, lend = 1, ljoin = 1)
rect(xleft = 0.2, ybottom = c(-0.3, 0.3), xright = 5, ytop = c(-1.2, 1.2),
     col = "white")

text(x = 1.1, y = 0.7, labels = "An author's work", pos = 4, adj = c(1, 1))
text(x = 0.9, y = -0.8, labels = "A reader's memory", pos = 4, adj=c(1, 1))
rect(xleft = c(0, 3, 6, 0, 3, 6), 
     ybottom = c(5.2, 3.4, 1.6, -5.2, -3.4, -1.6),  
     xright = c(3.5, 6.5, 9.5, 3.5, 6.5, 9.5), 
     ytop = c(6.6, 4.8, 3.0, -6.6, -4.8, -3.0), col = "white")
rect(xleft = 9.0, ybottom = -1.0, xright = 10.50, ytop = 1.20,
     col = "gray90", border = NULL)
rect(xleft = 8.9, ybottom = -1.1, xright = 10.35, ytop = 1.05, col="white")

text(x = c(0.1, 3.1, 6.1), y = c(5.8, 3.95, 2.15, -6.05, -4.25, -2.4), 
     pos = 4, adj = c(1, 1), labels = c(
       "1. Research Idea\n- One sentence\n- A few days or months",
       "2. Outline\n- 2 pages\n- Several months",
       "3. Details\n- 30 pages\n- Several weeks",
       "3. Research Idea\n- Topic only\n- After several years",
       "2. Outline\n- Structure only\n- In a year",
       "1. Details\n- Most facts\n- In several days"))
text(x = 9.1, y = 0.4, labels = "Paper", pos = 4, adj = c(1, 1))
memo.graph <- recordPlot()  # record the graph for version C


# aDemand -----------------------------------------------------------------

# 1. Create a screen device
win.graph(width = 5, height = 4, pointsize = 9); bringToTop(stay = TRUE)

# 2a. Parameter choices for version A
# in.mai <- c(0.5, 0.5, 0.1, 0.1)
# in.axes <- in.ann <- in.draw <- TRUE; in.x <- 1.5; in.y <- 1.5

# 2b. Parameter choices for version B
in.mai <- c(0.3, 0.3, 0.1, 0.1)
in.axes <- in.ann <- in.draw <- FALSE; in.x <- 0.2; in.y <- 0.3

# 3. Shared commands
par(mai = in.mai, mgp = c(2, 0.6, 0), las = 1, family = "serif")
plot(0, xlim = c(0, 10), ylim = c(1, 9.5), type = "n", axes = in.axes,
     ann = in.ann, xaxs = "i", yaxs = "i", xlab = "", ylab = "")

xa <- 4.5; xb <- 5.25; xc <- 4.5; ya <- 5.0; yb <- 5.50; yc <- 6.0
polygon(x = c(0, xa, xb, 0), y = c(ya, ya, yb, yb), border = NA,
        density = 30, lty = "dotted")  # producer surplus
polygon(x = c(0, xb, xc, 0), y = c(yb, yb, yc, yc), border = NA,
        col = "grey90")  # consumer surplus

if (in.draw) {  # adjustment on the axis lines
  box()
} else {
  axis(side = 1, labels = FALSE, lwd.ticks = -1, at = c(0, 10))
  axis(side = 2, labels = FALSE, lwd.ticks = -1, at = c(1, 9.5))
}
lines(x = c(0, 9), y = c(2, 8))  # Supply curve
lines(x = c(0, 9), y = c(8, 2))  # Initial demand curve
lines(x = c(0, 9), y = c(9, 3))  # New demand curve

lines(x = c(0,  xa), y = c(ya, ya), lty = 2)  # five dashed lines
lines(x = c(0,  xb), y = c(yb, yb), lty = 2)
lines(x = c(0,  xc), y = c(yc, yc), lty = 2)
lines(x = c(xc, xc), y = c(0,  yc), lty = 2)
lines(x = c(xb, xb), y = c(0,  yb), lty = 2)

points(x = c(xa, xb), y = c(ya, yb), pch = 16, cex = 1.5)  # 3 points
points(x = xc, y = yc, pch = 21, bg = "grey70")
text(x = c(xa - 0.6, xb, xc), y = c(ya + 0.1, yb, yc),
     labels = c("a", "b", "c"), pos = c(1,3,3))  # labels at plot region
text(x = c(9, 9, 9), y = c(2, 3, 8) + 0.2,
     labels = c(ex pression(D[a]), ex pression(D[b]), "S"), pos = 3 )

mtext(text = c(ex pression(Q[a]), ex pression(Q[b])),
      side = 1, line = in.x, at = c(xa, xb))  # labels in figure margin
mtext(text = c("d", "e", "f", ex pression(P[b]), ex pression(P[a]), "g"),
      side = 2, line = in.y,  at = c(9, 8, yc, yb, ya, 2))
out.demand <- recordPlot()

# Graph version C

Bibliography

citation(“erer”)

To cite package ‘erer’ in publications use:

Changyou Sun (2020). erer: Empirical Research in Economics with R. R package version 3.0.
https://CRAN.R-project.org/package=erer

LaTeX的用户的BibTeX条目是

@Manual{,
title = {erer: Empirical Research in Economics with R},
author = {Changyou Sun},
year = {2020},
note = {R package version 3.0},
url = {https://CRAN.R-project.org/package=erer},
}

ATTENTION: This citation information has been auto-generated from the package DEsc riptION file and may need
manual editing, see ‘help(“citation”)’.

扫码加我 拉你入群

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

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

分享0 收藏3 回帖

关键词：计量经济学 计量经济 经济学家 ERER CHAN

tulipsliu 发表于 2021-7-15 18:23
欢迎使用Markdown编辑器

经管之家：Do the best economic and management education！

好的好的好的好的

yangming98 发表于 2021-7-15 22:36
好的好的好的好的

哦？ 你的拼音像我以前认识的一个人的名字。
代码有用就好。

$$
\partial_s f(x) = \frac{\partial}{\partial x_0} f(x)\quad
\text{for $x= x_0 + I x_1$.}

$$

$$
\begin{align}
  \partial C&\!: &D_1U(C,1-L)&=\lambda\cdot(1-2\gamma V)\label{eq:d0C} \\
  \partial L&\!: &D_2U(C,1-L)&=-w\dot \mu -\dot w\mu
   +\left( {\beta +\delta } \right)\mu w \label{eq:d0L} \\
  \partial Y&\!: &{\lambda  \over P}&=-\dot \mu +\left( {\beta +\delta } \right)\mu
     \label{eq:d0Y} \\
  \partial B_C&\!: & -{{\dot \lambda } \over P}
    &=\left( {r-\beta -{{\dot P} \over P}} \right){\lambda  \over P} \label{eq:d0B}\\
  \partial M&\!: & -{{\dot \lambda } \over P}+\left( {\beta +{{\dot P} \over P}} \right)
       {\lambda  \over P}&=\gamma {\lambda  \over P}V^2 \label{eq:d0M}
\end{align}
$$

$$
\color{red}{\mathop{\arg\max\limits}_{\{\pi ,\ell ,L,Y,B_F\in\mathscr{B}\}}\int\limits_0^\infty
      {e^{-\beta t}\phi \left( {\pi _t} \right)dt}
}
$$

$$
\begin{align}
   \zeta&\!: &\pi &\le f(L)-{Y \over P}+{{\dot B_F} \over P}-{{rB_F} \over P}
   \label{eq:zeta0}\\
   \nu&\!: & \dot Y-w\dot L&\ge -\delta \left( {Y-wL} \right) \label{eq:nu0} \\
   && B &\geq 0 \:. \label{eq:B0ge0}
\end{align}
$$

$$
\begin{align}
   \mu (t)&=\int\limits_0^\infty  {e^{-\left( {\beta +\delta } \right)s}{{\lambda
   _{t+s}} \over {P_{t+s}}}ds} \label{eq:0muSol}\\
   \nu (t)&=\int\limits_0^\infty  {e^{-\left( {\beta +\delta } \right)s}{{\zeta _{t+s}}
   \over {P_{t+s}}}ds} \label{eq:0nuSol}
\end{align}
$$

$$
   {{D_2U\cdot \left( {1+2\gamma V} \right)} \over {D_1U}}={{-\dot \mu +\left( {\beta
   +\delta } \right)\mu } \over \lambda }
$$

$$
f'=-{{\dot w} \over w}{\nu  \over \zeta }+{{-\dot \nu w+\left( {\beta +\delta }
   \right)\nu w} \over \zeta }=-{{\dot w} \over w}{\nu  \over \zeta }+{w \over P}\:.
$$