[关闭]
@fanxy 2020-03-12T18:34:34.000000Z 字数 5439 阅读 3957

一、静态最优数学基础

樊潇彦 复旦大学经济学院 数量经济学


准备工作

  1. setwd("D:\\...") # 设定工作目录
  2. rm(list=ls())
  3. # 安装程序包
  4. install.packages("Rsolnp") # 非线性最优
  5. install.packages("ggplot2") # 作图
  1. # 调用
  2. library(Rsolnp)
  3. library(ggplot2)

1. 微积分

1.1 最优投资组合

  1. # 作图
  2. mu=0.1; r_f = 0.02; A =10; sigma=0.2
  3. theta= seq(0,1, by=0.05)
  4. U = theta*mu + (1-theta)*r_f - A*theta^2*sigma^2/2
  5. data=data.frame(theta, U)
  6. ggplot(data) +
  7. geom_line(aes(x=theta, y=U)) +
  8. geom_vline(xintercept=(mu-r_f)/A/sigma^2, color="red") +
  9. labs(x="风险资产占比", y="效用") +
  10. theme_bw()
  11. # 求解
  12. theta0=0.01; LB=0; UB=1 # 参数初始值、下限和上限
  13. U_fun=function(theta){ # 目标函数
  14. U=theta*mu + (1-theta)*r_f - A*theta^2*sigma^2/2
  15. return(-U) # 默认为最小化,加负号改为最大化
  16. }
  17. -U_fun(0.2) # 当theta=0.2时,计算效用值
  18. UMP=solnp(theta0, fun=U_fun, LB=LB, UB=UB) # 求解消费者最优化问题
  19. result=data.frame(theta=UMP$pars,maxU=-tail(UMP$values,1),SOC=UMP$hessian)
  20. round(result,3) # 保留小数点后3位,显示最优化结果

1.2 消费者最优

  1. ## 第1步:函数设定
  2. U_fun=function(c){
  3. if (rho==0) {U=(c[1])^alpha*(c[2])^(1-alpha)}else{
  4. U=(alpha*(c[1])^rho+(1-alpha)*(c[2])^rho)^(1/rho) # 最大化效用
  5. }
  6. return(-U) # min_(-U) # 默认为最小化目标函数
  7. }
  8. budget=function(c){p1*c[1]+p2*c[2]} # 预算约束方程
  9. c2_u0=function(c2,alpha,rho,c1,maxU0){ # 给定U0和c1求c2,用于画无差异曲线
  10. if (rho==0){log(maxU0)-( alpha*log(c1)+(1-alpha)*log(c2) )}else log(maxU0)-1/rho*log(alpha*c1^rho+(1-alpha)*c2^rho)
  11. }
  12. ## 第2步:分不同情况求解
  13. fig_data=result=data.frame() # 准备存放作图的数据
  14. for (case in c("Basic","Love_Apple","Expensive_Apple","High_W","CES")){ # 分四种情况
  15. p1=0.1; p2=1; w=300; # 外生变量
  16. alpha=0.3; sigma=1; rho=1-1/sigma # 模型参数
  17. switch(case,
  18. "Love_Apple"={
  19. alpha=0.5
  20. },
  21. "Expensive_Apple"={
  22. p1=0.15
  23. },
  24. "High_W"={
  25. w=330
  26. },
  27. "CES"={
  28. sigma=2; rho=1-1/sigma
  29. })
  30. maxc1=w/p1; maxc2=w/p2
  31. # 给定预算约束,求解消费者最优化问题
  32. find_cstar=solnp(c(maxc1/2,maxc2/2), # 初始猜测
  33. fun=U_fun, # 目标函数
  34. eqfun=budget, eqB=w, # 预算约束方程和收入
  35. LB=c(0.001,0.001), UB=c(maxc1,maxc2)) # 搜索下界和上界
  36. cstar=find_cstar$pars # 最优化结果
  37. maxU0=-U_fun(cstar) # 实际的最大化效用,用于后面画无差异曲线
  38. result=rbind(result,data.frame(case=case,c1=cstar[1],c2=cstar[2],maxU0=maxU0))
  39. c1vector=seq(0.1*maxc1,0.7*maxc1,length.out=10);
  40. c2vector=data.frame()
  41. for (c1 in c1vector){
  42. c2= uniroot(c2_u0, # 求无差异曲线上的c2
  43. c(0.001,1000), # 求解范围
  44. alpha=alpha,rho=rho,c1=c1,maxU0=maxU0,tol=0.0001)
  45. c2vector=rbind(c2vector,c2$root)
  46. }
  47. data=data.frame(case=case, U0=maxU0, c1=c1vector, c2=c2vector,
  48. c2_budget=(w-p1*c1vector)/p2) # 加上预期约束方程
  49. colnames(data)=c("case","U0","c1","c2","c2_budget") # 指标名称
  50. if (case=="Basic"){
  51. Basicdata=data[,3:5]
  52. colnames(Basicdata)=c("B_c1","B_c2","B_c2_budget") # 基准情况
  53. }else{
  54. fig_data=rbind(fig_data,cbind(Basicdata,data))
  55. }
  56. }
  57. # 报告最优化结果
  58. print(result)
  59. # 作图
  60. ggplot(fig_data,aes(c1,c2,color="red"))+geom_line()+geom_line(aes(c1,c2_budget,color="red"))+
  61. geom_line(aes(B_c1,B_c2,color="blue"))+geom_line(aes(B_c1,B_c2_budget,color="blue"))+
  62. facet_wrap(~case,scales="free")+ xlab("Apple")+ylab("Fish")+
  63. guides(color = guide_legend(title = NULL)) + theme_bw() +
  64. theme(legend.position = "non")

image_1ca2qc7vpu6i1d761g0gnoi1ukj9.png-47.9kB

2. 线性代数

2.1 矩阵运算的基本命令

  1. A=matrix(c(3,1,2,4),nrow=2) # 建立矩阵
  2. rowSums(A) # 行加总
  3. colSums(A) # 列加总
  4. a=A[,1]; t(a)%*%a # 第1列向量的欧氏范数
  5. I=diag(2) # 2维单位阵
  6. A_T = t(A) # 转置
  7. B=A+I # 加法
  8. solve(B) # 求逆
  9. A%*%B # 矩阵乘法

对应于矩阵的特征值 的列特征向量记为 为特征矩阵(以特征值为主对角线),则:


对应于矩阵的特征值 的行特征向量记为 ,相应地有:

由于 ,因此 的列特征向量矩阵的转置矩阵。

  1. lambda=eigen(A)$value # 特征值
  2. colP=eigen(A)$vectors # 特征向量,默认为列向量
  3. A%*%colP # 检验 AP=PB
  4. colP%*%diag(lambda)
  5. rowP=t(eigen(t(A))$vectors) # 特征行向量, 等于t(A)的特征列向量的转置
  6. rowP%*%A # 检验 rowP*A=B*rowP
  7. diag(lambda)%*%rowP

2.2 剑桥食谱

  1. A=matrix(c(36,51,13,
  2. 52,34,74,
  3. 0,7,1.1),nrow=3,byrow=T)
  4. b=c(33,45,3)
  5. solve(A,b)

2.3 两部门宏观经济模型

  1. c=0.9; beta=0.5
  2. K=matrix(c(1-c, -c,
  3. beta*(1-c), 1-beta*c),nrow=2,byrow=T)
  4. C0=100; I0=1000
  5. v=c(C0,I0+beta*C0)
  6. solve(K)%*%v

2.4 投入-产出表

  1. M=matrix(c(15,20,30,
  2. 30,10,45,
  3. 20,60,0),nrow=3,byrow=T) # 中间投入矩阵
  4. d=c(35,115,70) # 最终需求
  5. x=c(100,200,150) # 总产出
  6. A=M/matrix(rep(x,each=3),nrow=3) # 直接消耗系数矩阵
  7. solve(diag(3)-A)%*%d # 检验最终产出公式
  8. solve(diag(3)-A)%*%c(100,200,300) # 给定新的最终需求,求总产出

3. 统计与计量

3.1 《定价未来》

  • 1834年生于贝当古的法国人朱尔斯 雷格纳特(Jules Regnault)是我们故事中的一个被遗忘的英雄。...... 他的父亲是一位海关关员,在他12岁时就在穷困潦倒中去世,接着他家搬到了布鲁塞尔。作为穷人家的孩子,19岁的朱尔斯被免除了兵役,以当抄写员的收入来维持家用。与此同时,他旁听了布鲁塞尔自由大学的高级数学课程,同样被免除了学费,但他并没有拿到毕业证。

  • 当朱尔斯28岁时,他和他的弟弟奥迪朗搬到了巴黎。他俩合租了一个很小的公寓,这个公寓面积很小并且处于顶楼,因此他们不需缴税。他们的陋室还有个优点,就是非常靠近巴黎交易所。巴黎交易所成立于1862年,两兄弟都在交易所里给经纪人助理。当他们到来时,市场已经经过了一些年的平稳增长,即将进入一个高波动的时期,这将为投资者带来风险的同时创造机遇。

  • 当时的不稳定局面引起了金融界、法律界和经济学界的专家以及政府官员的争论,焦点集中于如何才能避免市场波动,实现可持续增长并回避灾难。......如何才能阻止投机者对金融市场功能有序发挥的破坏、促进经济持续增长呢?雷格纳特正是为了回答这个问题而完成了他的专著。

http://pierre.driout.perso.sfr.fr/Pages-148-149.jpg

3.2 用 lm() 做OLS回归

我们假定真实参数 ,并模拟出

  1. # 生成模拟数据
  2. N=200 # 样本数
  3. beta=c(1, 0.5, 0.1, 1) # 参数
  4. K=length(beta)
  5. x=matrix(rep(NaN,N*K),nrow=N) # N*K 空矩阵
  6. mu_x=c(1,0,0,0); sd_x=c(0,1,2,3) # 均值和标准差
  7. set.seed(1) # 设定随机种子使模拟结果可复制
  8. for (k in 1:K) { # 模拟生成x和y
  9. x[,k]=rnorm(N,mu_x[k],sd_x[k])
  10. }
  11. y=x%*%beta
  12. colnames(x)=c("con","x1","x2","epsilon")
  13. mydata=data.frame(y, x)
  14. # 估计
  15. x=x[,1:3]
  16. beta_bar=solve(t(x)%*%x)%*%t(x)%*%y # 线性最小二乘估计量
  17. myOLS=lm(y~x1+x2, data=mydata) # OLS 回归
  18. summary(myOLS) # 报告OLS回归结果

如果某个解释变量是离散变量(如年龄组),想把它设为哑变量,就用factor()。此外Kabacoff(2013)列出了 lm() 的其他常用设定,可以用来在解释变量中加入二次项、交互项,以及去掉某些解释变量:
lm.jpg-183.8kB

4. 数据与软件

4.1 数据来源

4.2 软件资料

添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注