EF_BookCode/Ch17_DecisonSupport.Rmd at master · EcoForecast/EF_BookCode · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
---
title: "Ch16 Decision Support"
author: "Michael Dietze"
date: "September 16, 2014"
output: html_document
---

```{r}

obs = c(9,7,4,3,1,1.25,1.5,1.75,2,-1,-1,-3,-6)
adjust=1.25
plot(density(obs,adjust=adjust))


filled.density <- function(x,adjust,vertical=TRUE,offset=0,...){
  den = density(x,adjust=adjust)
  if(vertical){
    polygon(c(0,den$y,0)+offset,c(den$x[1],den$x,den$x[length(den$x)]),...)
  }else{
    polygon(c(den$x[1],den$x,den$x[length(x)]),c(0,den$y,0)+offset,...)
  }
}

col.alpha <- function(col,alpha=1){
  rgb = col2rgb(col)
  rgb(rgb[1],rgb[2],rgb[3],alpha*255,maxColorValue=255)
}

##### define parameters
xlim = c(-8,10)
x = seq(xlim[1],xlim[2],length=1000)
ymax = max(density(obs,adjust=adjust)$y)
ymin = -ymax*0.2
ylim=c(ymin,ymax)

alpha=0.4
plot(density(obs,adjust=adjust),xlim=xlim,ylim=ylim,ylab="Density",type='n',bty='n',
       xaxt="n",yaxt="n",xlab="Land Uptake (GtC/yr)",mgp=c(1,1,0)*1.7,
     main=" ")
axis(1, pos=ymin)
axis(2, pos=x[1])
filled.density(obs,adjust,FALSE,col=col.alpha(4,alpha),lwd=2)

## decisions
dbreaks = c(xlim[1],-3,0,2,5,xlim[2])
cols = c(2,3,5,6,7)
for(i in 1:(length(dbreaks)-1)){
  polygon(c(dbreaks[i],dbreaks[i+1],dbreaks[i+1],dbreaks[i]),
          c(ymin,ymin,0,0),col=col.alpha(cols[i],0.4))
  text((dbreaks[i]+dbreaks[i+1])/2,ymin/2,i,cex=1.5)
}

## Bringing in risk tolerance explicitly
col.a = tapply(cols,seq_along(cols),col.alpha,0.4)
n = 18
risk = seq(-1,1,length=7)
sink = seq(-8,10,length=n)
dec.matrix = matrix(c(2,2,2,2,3,3,4,4,4,5,5,5,5,5,5,5,5,5,
                      1,1,2,2,2,3,3,4,4,4,4,5,5,5,5,5,5,5,
                      1,1,1,1,2,2,2,3,3,4,4,4,5,5,5,5,5,5,
                      1,1,1,1,1,2,2,2,3,3,4,4,4,5,5,5,5,5,
                      1,1,1,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,
                      1,1,1,1,1,1,1,1,2,2,2,3,3,3,4,4,4,4,
                      1,1,1,1,1,1,1,1,1,1,2,2,2,3,3,3,3,4
                    ),7,n,byrow=TRUE)
image(sink,risk,t(dec.matrix),col=col.a,xlab="Land Uptake (GtC/yr)",ylab="Risk Adversion",cex.axis=1.3,cex.lab=1.5)

```


## UTILITY

```{r}
y = 10^seq(1.5,3.2,length=1000)
Y = c(300,1000)
U = c(1,2)
m = diff(U)/diff(log(Y))
b = U[1]-m*log(Y[1])
u = m*log(y)+b
plot(y,u,xlab="Dollars ($)",ylab="Utility",
     type='l',lwd=3,ylim=c(0,2.5),xlim=c(0,1500),
     cex.lab=1.4,cex.axis=1.2)
points(Y,U,cex=2)
lines(c(-100,Y[1],Y[1]),c(U[1],U[1],-3),lty=3,lwd=3)
lines(c(-100,Y[2],Y[2]),c(U[2],U[2],-3),lty=3,lwd=3)

```