索洛模型的R語言簡單實現

01-24

中宏課的一次作業，老師讓用Matlab實現索洛模型，但我之前沒怎麼接觸過Matlab，就用R語言簡單得寫了一個，分享下。

模型

索洛模型的關鍵方程

$dot{k_t}=sf(k_t)-(n+g+delta)k_t$

生產函數C-D型的關鍵方程

$dot{k_t}=sk_t^alpha-(n+g+delta)k_t$

函數定義與求解

```{r plot-1, echo=TRUE}library(deSolve)library(ggplot2)# 定義函數fun <- function(Time, State, Pars) { with(as.list(c(State, Pars)), { dk <- s*(k^a)-(n+g+d)*k return(list(c(dk))) })}# 賦參數值pars <- c(s = 0.6, n=0.02, g=0.05, d=0.5, a=0.5)# 時間階段（下標t）times <- seq(0, 50, by = 1)# 設定初值k_ini_1 <- c(k = 0.01)k_ini_2 <- c(k = 3)# 常微分方程求解out1 <- as.data.frame(ode(k_ini_1, times, fun, pars))out1$cl <- "1"out2 <- as.data.frame(ode(k_ini_2, times, fun, pars))out2$cl <- "2"out <- rbind(out1, out2)out_display <- cbind(out1, out2)[,c(1,2,4,5)]print(out_display)```

k隨迭代次數變化趨勢

```{r}# 圖 1 k隨迭代次數變化趨勢ggplot(data = out, aes(x = time , y = k, shape=cl, lty=cl, col=cl)) + geom_point() + geom_line() + scale_shape_discrete(breaks=c("1", "2"), labels=c("k初始值為0.01", "k初始值為3.00"))+ scale_color_discrete(breaks=c("1", "2"), labels=c("k初始值為0.01", "k初始值為3.00"))+ scale_linetype_discrete(breaks=c("1", "2"), labels=c("k初始值為0.01", "k初始值為3.00"))+ labs(title="k隨迭代次數變化趨勢")+ theme(legend.title = element_blank())```

實際投資與持平投資

```{r}# 圖 2 實際投資與持平投資library(reshape2)out$sfk <- 0.6 * (out$k)^0.5out$ngdk <- (0.02 + 0.05 + 0.5) * out$kout_melt <- melt(data = out, id.vars = c("time", "k", "cl"))ggplot(data = out_melt, aes(x = k , y = value, shape=variable, lty=variable, col=variable)) + geom_point() + geom_line() + scale_shape_discrete(breaks=c("sfk", "ngdk"), labels=c("實際投資", "持平投資"))+ scale_color_discrete(breaks=c("sfk", "ngdk"), labels=c("實際投資", "持平投資"))+ scale_linetype_discrete(breaks=c("sfk", "ngdk"), labels=c("實際投資", "持平投資"))+ labs(title="實際投資與持平投資",y="每單位有效勞動的平均投資")+ theme(legend.title = element_blank())```

索洛模型中k的相圖

```{r}# 圖 3 索洛模型中k的相圖out$kdot <- 0.6 * out$k^0.5 - (0.02 + 0.05 + 0.5) * out$kggplot(data = out, aes(x = k , y = kdot, shape=cl, lty=cl, col=cl)) + geom_point() + geom_line() + geom_hline(yintercept = 0) + scale_shape_discrete(breaks=c("1", "2"), labels=c("k初始值為0.01", "k初始值為3.00"))+ scale_color_discrete(breaks=c("1", "2"), labels=c("k初始值為0.01", "k初始值為3.00"))+ scale_linetype_discrete(breaks=c("1", "2"), labels=c("k初始值為0.01", "k初始值為3.00"))+ labs(title="索洛模型中k的相圖")+ theme(legend.title = element_blank())```