@fanxy
2016-07-08T15:57:30.000000Z
字数 3162
阅读 2446
樊潇彦
复旦大学经济学院
经济数学
# 准备工作
setwd("D:\\...")
install.packages("Rdonlp2", repos="http://R-Forge.R-project.org")
install.packages("lpSolve")
install.packages(c("ggplot2","plotly"))
library(lpSolve)
library(Rdonlp2)
library(ggplot2)
library(plotly)
# 生产计划
production = lp(objective.in=c(24, 18),
const.mat=matrix(c(1, 1, 2, 3, 3, 2, 1, 0, 0,1), nrow=5, byrow=T),
const.rhs=c(150, 240, 300,0,0),
const.dir=c("<=", "<=", "<=",">=",">="),
direction="max")
production$solution # 规划结果
production$objval # 目标值
# 运输问题
transportation = lp(objective.in=c(4,5,6,5,2,4),
const.mat=rbind(matrix(c(1, 1, 1, 0, 0, 0,
0, 0, 0, 1, 1, 1,
1, 0, 0, 1, 0, 0,
0, 1, 0, 0, 1, 0,
0, 0, 1, 0, 0, 1), nrow=5,byrow=T),
diag(6)),
const.rhs=c(12, 8, 8,6,6,rep(0,6)),
const.dir=c("<=", "<=", ">=",">=",">=",
">=",">=",">=", ">=",">=",">="),
direction="min")
transportation$solution # 规划结果
transportation$objval # 目标值
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()
r=0.1; w=5; alpha=0.5; Cost=100; N=12
K=seq(0.01,Cost/r, length.out=N)
L=seq(0.01,Cost/w, length.out=N)
# 用ggplot2包做二维图
firmdata <- transform(expand.grid(K=K,L=L), Y=(K*L)^alpha)
ggplot(firmdata) +
stat_contour(aes(x=K, y=L, z=Y, colour=Y)) +
geom_abline(slope=-r/w, intercept=Cost/w, color="red") +
labs(xlab="K", ylab="L") +
theme_bw()
# 用plotly包作三维图
Y=firmdata$Y
dim(Y)=c(N,N)
plot_ly(z=Y, type="contour") %>% # 两维产出等高线
layout(scene = list(
xaxis = list(title = "K"),
yaxis = list(title = "L"),
zaxis = list(title = "Y")))
plot_ly(z=Y, type="surface") %>% # 三维产出曲面
layout(scene = list(
xaxis = list(title = "K"),
yaxis = list(title = "L"),
zaxis = list(title = "Y")))
# 求解
par.l= c(0,0); par.u = c(100/r,100/w) # 参数的下限和上限
fn=function(x){ # 目标函数
f=(x[1]*x[2])^alpha-r*x[1]-w*x[2]
return(-f) # 默认为最小化,加负号改为最大化
}
A=matrix(c(r,w),nrow=1) # 线性约束
lin.l=c(Cost); lin.u=c(Cost) # 线性约束上限和下限
p0=c(1,1) # 参数(K和L)初始值
ret=donlp2(p0,fn,par.l=par.l,par.u=par.u,A,lin.l=lin.l,lin.u=lin.u) # 解最优化问题
K=ret$par[1]; L=ret$par[2]
data.frame(K,L) # 最优化结果
alpha=0.4; p1=10; p2=5; inc=100; N=12
# 作图
c1=seq(0.01, inc/p1, length.out=N)
c2=seq(0.01, inc/p2, length.out=N)
data <- transform(expand.grid(c1=c1,c2=c2),
U=alpha*log(c1)+(1-alpha)*log(c2))
U=data$U
dim(U)=c(N,N)
plot_ly(z=U, type="contour") %>%
layout(scene = list(
xaxis = list(title = "C1"),
yaxis = list(title = "C2"),
zaxis = list(title = "U")))
# 求解
par.l= c(0,0); par.u = c(inc/p1,inc/p2)
u=function(x){
u=alpha*log(x[1])+(1-alpha)*log(x[2])
return(-u)
}
A=matrix(c(p1,p2),1)
lin.l=c(inc); lin.u=c(inc)
p0=c(1,1)
ret=donlp2(p0,fn=u,par.l=par.l,par.u=par.u,A,lin.l=lin.l,lin.u=lin.u)
c1=ret$par[1]; c2=ret$par[2]
data.frame(c1,c2)
# 结果检查
c1_check=alpha*inc/p1
c2_check=(1-alpha)*inc/p2
data.frame(c1_check,c2_check)