@fanxy
2020-03-12T18:34:34.000000Z
字数 5439
阅读 3957
樊潇彦
复旦大学经济学院
数量经济学
setwd("D:\\...") # 设定工作目录
rm(list=ls())
# 安装程序包
install.packages("Rsolnp") # 非线性最优
install.packages("ggplot2") # 作图
# 调用
library(Rsolnp)
library(ggplot2)
# 作图
mu=0.1; r_f = 0.02; A =10; sigma=0.2
theta= seq(0,1, by=0.05)
U = theta*mu + (1-theta)*r_f - A*theta^2*sigma^2/2
data=data.frame(theta, U)
ggplot(data) +
geom_line(aes(x=theta, y=U)) +
geom_vline(xintercept=(mu-r_f)/A/sigma^2, color="red") +
labs(x="风险资产占比", y="效用") +
theme_bw()
# 求解
theta0=0.01; LB=0; UB=1 # 参数初始值、下限和上限
U_fun=function(theta){ # 目标函数
U=theta*mu + (1-theta)*r_f - A*theta^2*sigma^2/2
return(-U) # 默认为最小化,加负号改为最大化
}
-U_fun(0.2) # 当theta=0.2时,计算效用值
UMP=solnp(theta0, fun=U_fun, LB=LB, UB=UB) # 求解消费者最优化问题
result=data.frame(theta=UMP$pars,maxU=-tail(UMP$values,1),SOC=UMP$hessian)
round(result,3) # 保留小数点后3位,显示最优化结果
## 第1步:函数设定
U_fun=function(c){
if (rho==0) {U=(c[1])^alpha*(c[2])^(1-alpha)}else{
U=(alpha*(c[1])^rho+(1-alpha)*(c[2])^rho)^(1/rho) # 最大化效用
}
return(-U) # min_(-U) # 默认为最小化目标函数
}
budget=function(c){p1*c[1]+p2*c[2]} # 预算约束方程
c2_u0=function(c2,alpha,rho,c1,maxU0){ # 给定U0和c1求c2,用于画无差异曲线
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)
}
## 第2步:分不同情况求解
fig_data=result=data.frame() # 准备存放作图的数据
for (case in c("Basic","Love_Apple","Expensive_Apple","High_W","CES")){ # 分四种情况
p1=0.1; p2=1; w=300; # 外生变量
alpha=0.3; sigma=1; rho=1-1/sigma # 模型参数
switch(case,
"Love_Apple"={
alpha=0.5
},
"Expensive_Apple"={
p1=0.15
},
"High_W"={
w=330
},
"CES"={
sigma=2; rho=1-1/sigma
})
maxc1=w/p1; maxc2=w/p2
# 给定预算约束,求解消费者最优化问题
find_cstar=solnp(c(maxc1/2,maxc2/2), # 初始猜测
fun=U_fun, # 目标函数
eqfun=budget, eqB=w, # 预算约束方程和收入
LB=c(0.001,0.001), UB=c(maxc1,maxc2)) # 搜索下界和上界
cstar=find_cstar$pars # 最优化结果
maxU0=-U_fun(cstar) # 实际的最大化效用,用于后面画无差异曲线
result=rbind(result,data.frame(case=case,c1=cstar[1],c2=cstar[2],maxU0=maxU0))
c1vector=seq(0.1*maxc1,0.7*maxc1,length.out=10);
c2vector=data.frame()
for (c1 in c1vector){
c2= uniroot(c2_u0, # 求无差异曲线上的c2
c(0.001,1000), # 求解范围
alpha=alpha,rho=rho,c1=c1,maxU0=maxU0,tol=0.0001)
c2vector=rbind(c2vector,c2$root)
}
data=data.frame(case=case, U0=maxU0, c1=c1vector, c2=c2vector,
c2_budget=(w-p1*c1vector)/p2) # 加上预期约束方程
colnames(data)=c("case","U0","c1","c2","c2_budget") # 指标名称
if (case=="Basic"){
Basicdata=data[,3:5]
colnames(Basicdata)=c("B_c1","B_c2","B_c2_budget") # 基准情况
}else{
fig_data=rbind(fig_data,cbind(Basicdata,data))
}
}
# 报告最优化结果
print(result)
# 作图
ggplot(fig_data,aes(c1,c2,color="red"))+geom_line()+geom_line(aes(c1,c2_budget,color="red"))+
geom_line(aes(B_c1,B_c2,color="blue"))+geom_line(aes(B_c1,B_c2_budget,color="blue"))+
facet_wrap(~case,scales="free")+ xlab("Apple")+ylab("Fish")+
guides(color = guide_legend(title = NULL)) + theme_bw() +
theme(legend.position = "non")
A=matrix(c(3,1,2,4),nrow=2) # 建立矩阵
rowSums(A) # 行加总
colSums(A) # 列加总
a=A[,1]; t(a)%*%a # 第1列向量的欧氏范数
I=diag(2) # 2维单位阵
A_T = t(A) # 转置
B=A+I # 加法
solve(B) # 求逆
A%*%B # 矩阵乘法
对应于矩阵的特征值 的列特征向量记为 , 为特征矩阵(以特征值为主对角线),则:
lambda=eigen(A)$value # 特征值
colP=eigen(A)$vectors # 特征向量,默认为列向量
A%*%colP # 检验 AP=PB
colP%*%diag(lambda)
rowP=t(eigen(t(A))$vectors) # 特征行向量, 等于t(A)的特征列向量的转置
rowP%*%A # 检验 rowP*A=B*rowP
diag(lambda)%*%rowP
A=matrix(c(36,51,13,
52,34,74,
0,7,1.1),nrow=3,byrow=T)
b=c(33,45,3)
solve(A,b)
c=0.9; beta=0.5
K=matrix(c(1-c, -c,
beta*(1-c), 1-beta*c),nrow=2,byrow=T)
C0=100; I0=1000
v=c(C0,I0+beta*C0)
solve(K)%*%v
M=matrix(c(15,20,30,
30,10,45,
20,60,0),nrow=3,byrow=T) # 中间投入矩阵
d=c(35,115,70) # 最终需求
x=c(100,200,150) # 总产出
A=M/matrix(rep(x,each=3),nrow=3) # 直接消耗系数矩阵
solve(diag(3)-A)%*%d # 检验最终产出公式
solve(diag(3)-A)%*%c(100,200,300) # 给定新的最终需求,求总产出
1834年生于贝当古的法国人朱尔斯 雷格纳特(Jules Regnault)是我们故事中的一个被遗忘的英雄。...... 他的父亲是一位海关关员,在他12岁时就在穷困潦倒中去世,接着他家搬到了布鲁塞尔。作为穷人家的孩子,19岁的朱尔斯被免除了兵役,以当抄写员的收入来维持家用。与此同时,他旁听了布鲁塞尔自由大学的高级数学课程,同样被免除了学费,但他并没有拿到毕业证。
当朱尔斯28岁时,他和他的弟弟奥迪朗搬到了巴黎。他俩合租了一个很小的公寓,这个公寓面积很小并且处于顶楼,因此他们不需缴税。他们的陋室还有个优点,就是非常靠近巴黎交易所。巴黎交易所成立于1862年,两兄弟都在交易所里给经纪人助理。当他们到来时,市场已经经过了一些年的平稳增长,即将进入一个高波动的时期,这将为投资者带来风险的同时创造机遇。
当时的不稳定局面引起了金融界、法律界和经济学界的专家以及政府官员的争论,焦点集中于如何才能避免市场波动,实现可持续增长并回避灾难。......如何才能阻止投机者对金融市场功能有序发挥的破坏、促进经济持续增长呢?雷格纳特正是为了回答这个问题而完成了他的专著。
我们假定真实参数 ,并模拟出 和 :
# 生成模拟数据
N=200 # 样本数
beta=c(1, 0.5, 0.1, 1) # 参数
K=length(beta)
x=matrix(rep(NaN,N*K),nrow=N) # N*K 空矩阵
mu_x=c(1,0,0,0); sd_x=c(0,1,2,3) # 均值和标准差
set.seed(1) # 设定随机种子使模拟结果可复制
for (k in 1:K) { # 模拟生成x和y
x[,k]=rnorm(N,mu_x[k],sd_x[k])
}
y=x%*%beta
colnames(x)=c("con","x1","x2","epsilon")
mydata=data.frame(y, x)
# 估计
x=x[,1:3]
beta_bar=solve(t(x)%*%x)%*%t(x)%*%y # 线性最小二乘估计量
myOLS=lm(y~x1+x2, data=mydata) # OLS 回归
summary(myOLS) # 报告OLS回归结果
如果某个解释变量是离散变量(如年龄组),想把它设为哑变量,就用factor()
。此外Kabacoff(2013)列出了 lm()
的其他常用设定,可以用来在解释变量中加入二次项、交互项,以及去掉某些解释变量: