はじめに

このページでは、コブ＝ダグラス型の効用関数\(u=x^{\frac{1}{3}}y^{\frac{2}{3}}\)をRを用いて3Dで描画する方法を色々試します。

事前準備

パッケージの読み込み

今回用いるパッケージは以下の3つです。 rglとplotlyはプロットのインタラクティブな操作が可能です。

  • rgl
  • plot3D
  • plotly
library(rgl)
knitr::knit_hooks$set(webgl = hook_webgl)
library(plot3D)
library(plotly)
library(tidyverse)

効用関数の設定

ベクトルxとyで、各軸の値域とその刻み（今回は0.1刻み）を指定します。次に、xとyの組み合わせから直積を作った上で、それを先ほど定式化した効用関数に引数として与え、返り値を行列uとします。言い換えれば、xとyの組み合わせについて効用関数の値が計算された結果が行列uです。

x <- seq(0,6,by=0.1)
y <- seq(0,6,by=0.1)
u <- outer(x,y,function(x,y){(x^(1/3))*(y^(2/3))})

plot3Dを用いたプロット

persp3D関数を用い、上で作成したベクトル・行列を引数として効用関数をプロットします。プロットは、引数zで与えられた行列（今回は効用関数の値u）の要素の値に応じて色分けされます。

plot3D::persp3D(x=x,y=y,z=u)

rglを用いたプロット

persp3d関数を用い、先ほどとほぼ同じ設定でプロットできます。デフォルトの設定では、プロットの色は黒に描画されます。

rgl::persp3d(x=x,y=y,z=u)

色は引数colorで、不透明度はalphaで調整可能です。

rgl::persp3d(x=x,y=y,z=u,color="blue",alpha=0.4)

作成したプロットは、close3d関数で閉じます。

rgl::close3d()

plot3Dを用いた場合と同じような色分けをしたい場合には、パレットを手動で設定する必要があります。以下では、rainbowというパレットから、青〜赤に変化する100色を取ってきて描画用のパレットを作成しています。

color=rev(rainbow(100,start=0/6,end=4/6))
ucol=cut(u,100)

上のパレットを用いて効用関数をプロットします。

rgl::persp3d(x=x,y=y,z=u,color=color[ucol],alpha=0.8)

プロットを回転させるアニメーションを書き出すことも可能です。その際、事前にImageMagickというソフトをインストールしておく必要があります。

rgl::play3d(rgl::spin3d(axis=c(0,0,1)),duration=10)
rgl::movie3d(rgl::spin3d(axis=c(0,0,1)),duration=10,dir=getwd(),webshot=FALSE)

ミクロ経済学の教科書によくあるような、異なる効用水準で無差別曲線がどうなるかを説明するプロットを作ってみます。まず、surface3d関数を用いて、\(u=2\)と\(u=4\)の平面をプロットします。

rgl::surface3d(x=x,y=y,z=matrix(2,ncol=length(x),nrow=length(x)),col="black",alpha=0.8,add=T)
rgl::surface3d(x=x,y=y,z=matrix(4,ncol=length(x),nrow=length(x)),col="ivory",alpha=0.4,add=T)

\(u=4\)の場合は無差別曲線はx-y平面上に\(4=x^{\frac{1}{3}}y^{\frac{2}{3}} \rightarrow y = \left(\frac{4}{x^{1/3}} \right)^{3/2}\)、\(u=2\)の場合は\(2=x^{\frac{1}{3}}y^{\frac{2}{3}} \rightarrow y = \left(\frac{2}{x^{1/3}} \right)^{3/2}\)と描けるので、x軸の値域に対応したy軸方向の値を計算します。

u_4 <- cbind(seq(0.1,6,by=0.1),sapply(seq(0.1,6,by=0.1),function(x){(4/x^(1/3))^(3/2)}))
u_2 <- cbind(seq(0.1,6,by=0.1),sapply(seq(0.1,6,by=0.1),function(x){(2/x^(1/3))^(3/2)}))

それぞれの無差別曲線を2Dでプロットします。

plot(u_4,type="l",xlim=c(0,6),ylim=c(0,6),col="gray",xlab="x",ylab="y")
lines(u_2,type="l",col="black")
legend(4,6,legend=c("u=4","u=2"),col=c("gray","black"),lty=c(1,1),box.lty=0)

plotlyを用いたプロット

plot_ly関数を用い、これまでと同じ引数の設定でプロットできます。ただし、add_surface関数を同時に実行する必要があります。

fig <- plotly::plot_ly(x=x,y=y,z=u)
fig <- fig %>% 
  plotly::add_surface()
fig