---
title: "Workshop 1: Introduction to R Programming"
author: "oran.k@chula.ac.th"
date: "July 2022"
output: html_document
params:
  intab:  TRUE
---

```{r,setup,echo=F,include=FALSE}
  #rm(list=ls())
  knitr::opts_chunk$set(echo=TRUE,eval=TRUE,result="show"
                        ,warning = FALSE,message=FALSE)
  
  knitr::opts_template$set(hidden = list(
    fig.width = 10, fig.height = 6,
    fig.retina = 2, out.width = '100%',
    echo=!params$intab, eval=!params$intab,
    warning=FALSE,message=FALSE
  ))
  
  knitr::opts_template$set(badcode = list(
    echo=TRUE, eval=FALSE,
    warning=FALSE,message=FALSE
  ))
  
  #intab <- TRUE
```

```{r klippy, echo=FALSE, include=TRUE}
klippy::klippy(position = c('top', 'right'),color = 'darkred')
```

```{css, echo=F}
  .badcodeStyle { background-color: lightpink; font-family: courier }
  
```

```{css, echo=F}
  .hiddenStyle { background-color: lightblue; font-family: courier }
```

# Instruction `r if (params$intab) "{.tabset}"`

Read the [instruction](http://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/workshop01.pdf) carefully and think about how to develop R code to answer each questions.

------------------------------------------------------------------------

## Question 1 `r if (params$intab) "{.tabset}"`

##### Using R and RStudio to Monte Carlo Simulate: a summation of rolling dice experiment

### 01: create

##### create a distribution of rolling three fair-dices that has following description

-   All three dices are fair
-   The result of one replication of experiment is summation of three fair dices
-   There is 1000 replications

```{r  }
  threeDiceSum <- function(n){
    ## n<- 2
    return(ceiling(6*runif(n))+ceiling(6*runif(n))+ceiling(6*runif(n)))
  }

  set.seed(17)
  diceExp <- threeDiceSum(1000)
  hist(diceExp,col="grey",freq=F)
```

##### improve figure

```{r}
  hist(diceExp,col="grey",freq=F,main="Histrogram of Dice Experiment",
         breaks = seq(2.5,18.5,1),
         ylab="Freq", xlab="Sum of Values",
         cex =0.8)
  rug(jitter(diceExp),col="blue")
  axis(1,at=seq(3,18,1))
  box()
  
  abline(v=mean(diceExp),col="red")
  abline(v=mean(diceExp)+sd(diceExp),col="red",lty=2)
  abline(v=mean(diceExp)-sd(diceExp),col="red",lty=2)
  
  myText <- paste0("mean= ",format(mean(diceExp),digits = 3)," sd= ", format(sd(diceExp),digits = 3) )
  text(16,0.1,myText,cex=0.8)
  
  normApprox <- function(x){ dnorm(x,mean=mean(diceExp),sd=sd(diceExp) )}
  curve(normApprox,col="green",from=3,to=18,add=TRUE)
  
```

### 02: desp stat

##### what are the descriptive statistic?

-   **summary** min, mean, median, max, and sd

```{r}
  summary(diceExp)

  ##-- quantile/ 
  quantile(diceExp)
  
  ##-- percentile at any prob
  quantile(diceExp,prob=seq(0.1,0.9,0.1))
  
  ##-- frequency table
  table(diceExp)
```

### 03: stat test

##### Perform Shapiro Test (`shapiro.test()`) as well as Anderson Darling Test (`nortest::ad.test()`) of the previous experiment. The latter require 'nortest' package

{r} shapiro.test(diceExp)

nortest::ad.test(diceExp)

p-value is `5.337e-08` (which is \< 0.05) , we can not conclude that this is normal. Should we increase or decrease number of time experiment was performed

```{r hidden,echo=!params$intab,include=!params$intab }
  set.seed(17)

  diceExp <- threeDiceSum(100)
  shapiro.test(diceExp)
```

**Question**

-   Why test result better as number of replication reduce?

### EXTRA

##### create an unfair dice that:

-   a chance to roll any odd numbers are equal
-   a chance to roll any even numbers are equal
-   a chance to roll an odd number is twice as likely to an even number

```{r, opts.label="badcode",class.source="badcodeStyle"}
  ##-- This is incomplete code block. It is UNABLE TO EXECUTE 
  unFairDice <- function(n){ ??? }
  
  threeUnFairDiceSum <- function(n){
    ## n<- 2
    return(unFairDice(n)+unFairDice(n)+unFairDice(n) )
  }
```

```{r,opts.label="hidden",class.source="hiddenStyle"}
  unFairProb <- rep( c(2/9,1/9),3)
  unFairDice <- function(n){ sample(1:6,size = n,prob = unFairProb,replace=TRUE) }
  
  threeUnFairDiceSum <- function(n){
    ## n<- 2
    return(unFairDice(n)+unFairDice(n)+unFairDice(n) )
  }
  
  set.seed(17)
  result <- threeUnFairDiceSum(1000)
  hist(result,col="grey",freq=F)
  
  ##-- test for known function and AIC/BIC/LogLik
  broom::glance(MASS::fitdistr(result,"Poisson"))
  broom::glance(MASS::fitdistr(result,"lognormal"))
```

------------------------------------------------------------------------

## Question 2 \## `r if (params$intab) "{.tabset}"`

##### Using R and RStudio to create multiplication table $20 \times 20$

------------------------------------------------------------------------

### 01: loop

```{r}
### SOLUTION TO QUESTION 1 ###
tableA <- matrix(rep(NA,20*20),ncol=20) ##-- create dummy
for(i in 1:20){
  for(j in 1:20){
    tableA[i,j] <- i*j ##-- filled dummy with value
  }
}
head(tableA)
```

### 02: 'outer()'

##### Alternatively, we can use R function `outer()`

```{r}
tableB <- outer(1:20,1:20) ##-- use specific R cmd
head(tableB)
```

### EXTRA `r if (params$intab) "{.tabset}"`

#### tables

##### Can you do table of other functions, with `outer()`?

-   Additive table, i.e., $(a,b) = a+b$
-   Subtraction table, i.e., $(a,b) = a-b$
-   Power table, i.e, $(a,b) = a^b$

```{r,echo=!params$intab,include=!params$intab}
##-- set up function
addFn <- function(x,y){ x + y}
subFn <- function(x,y){ x - y}
powFn <- function(x,y){ x^y}

##-- use specific R cmd
tableC <- outer(1:20,1:20,powFn)
tableC[1:4,1:12]
```

#### access

##### Can access the matrix of multiplication table $20 \times 20$ in following condition?

-   access diagonal of matrix
-   access upper triangle of matrix
-   access reverse diagonal of matrix

```{r,echo=!params$intab,include=!params$intab}
##-- diagonal
result.1 <- tableA[diag(20) == 1]*diag(20)
result.1[1:4,1:12]


##-- upper triangle
result.2 <- tableA*matrix(as.numeric(upper.tri(tableA) ),ncol=20)
result.2[1:4,1:12]

##-- lower triangle 
result.2 <- tableA*matrix(as.numeric(lower.tri(tableA) ),ncol=20)


##-- reverse diagonal 
result.3 <- tableA[diag(20) == 0]*( (diag(20)+1) %% 2)
result.3[1:4,1:12]

```

## Question 3 \## `r if (params$intab) "{.tabset}"`

##### Determine **Net Present Value**, **Payback Period**, and **Internal Rate of Return** of the following cash flow of project if an annual investment interest rate is 8%

| Year     | 0     | 1   | 2   | 3   | 4   | 5   | 6    |
|----------|-------|-----|-----|-----|-----|-----|------|
| Cashflow | -1350 | 363 | 551 | 681 | 761 | 821 | 1467 |

------------------------------------------------------------------------

### 01: NPV

##### a) compute NPV

```{r}
### SOLUTION TO QUESTION 2A ###
cashflow <- c(-1350, 363, 551, 681, 761, 821, 1467)
intDef   <- 0.08

findPV <- function(cashflow, intDef){ ##-- function for PV
	period <- 0:(length(cashflow) -1)
	result <- cashflow/(1+intDef)^(period) 
	return(result)
}

findNPV <- function(cashflow, intDef){ ##-- function for NPV
  sum(findPV(cashflow, intDef)) 
} 

findNPV(cashflow, intDef) 

```

### 02: Payback

##### b) compute payback period

-   create cummerative present value

```{r}
### SOLUTION TO QUESTION 2B ###
periods <- 0:(length(cashflow) -1)  ##-- set periods

##-- calculate cumerative PV
cashflow.pv.sum <- cumsum(findPV(cashflow, intDef)) 
cashflow.pv.sum
```

-   show figure

```{r}
##-- plot figure----
plot(cashflow.pv.sum,periods,type="b",pch=16)
abline(v=0,col="red")
```

-   calculate value using R linear approximation function, called `approx()`

```{r}
##-- find interaction period when NPV = 0
payback <- approx(x=cashflow.pv.sum,y=periods,xout=0.0)  

payback$y
```

### 03: IRR

##### c) compute IRR

There are three methods to achieve the result:

##### Method 1: use polynomial by imaging

$$
0 = \frac{-1350}{(1+i)^0}+ \frac{363}{(1+i)^1} + \frac{551}{(1+i)^2} + \frac{681}{(1+i)^3} +  \frac{761}{(1+i)^4}+ \frac{821}{(1+i)^5} +\frac{1467}{(1+i)^6}
$$

Let $x^n = \frac{1}{(1+i)^n}$, the pervious equation becomes

$$
0 = -1350x^0+ 363x^1 + 551x^2 + 681x^3 +  761x^4+ 821x^5 +1467x^6
$$

we can find the root of $x$ by `polyroot()`

```{r}
root <- polyroot(cashflow)
root
```

the possible root should be positive and real, hence

```{r}
Re(root[4]) 
```

And, the solution is $i = \frac{1}{x} - 1$

```{r}
1/Re(root[4]) -1
```

##### Method 2: use figure and interpolation. this is similar with **payback period**

```{r}
### SOLUTION TO QUESTION 2C ###
interest   <- seq(0.0,1.0,0.01)
npv.sapply <- sapply(interest,findNPV,cashflow=cashflow)   ##-- use sapply to calculate 
plot(npv.sapply,interest,pch=16,cex=0.5,col="green")
abline(v=0,col="red") 

irr.appx <- approx(npv.sapply,interest,xout=0)
irr.appx$y
```

##### Method 3: use bi-section search function in R, called `uniroot()`.

```{r}
##-- find root (arg x : f(x) = 0) or IRR
irr.uni <- uniroot(findNPV,interval=c(0,1),cashflow=cashflow) 

irr.uni$root
```

##### **Discuss:** How to choose each method?

-   quality of solution
-   time to implement
-   limitation

------------------------------------------------------------------------

## Question 4 \## `r if (params$intab) "{.tabset}"`

##### A 10-year, 12% semiannual coupon bond, with a par value of 1000 USD may be called in 4 years at a call price of 1060 USD. The bond sells for 1100 USD (Assume that the bond has just been issued.)

------------------------------------------------------------------------

### 01: Maturity

##### a) What is the bond's yield to maturity?

```{r}
### SOLUTION TO QUESTION 3 ###
bondMat.cf<- c(-1100,rep(60,19),1060)
yield  <- uniroot(findNPV,interval=c(0,1),cashflow=bondMat.cf) 
yield$root
```

### 02: Call

##### b) What is the bond's yield to call?

```{r}
### SOLUTION TO QUESTION 3 ###
bondCall.cf<- c(-1100,rep(60,7),1060)
yield <- uniroot(findNPV,interval=c(0,1),cashflow=bondCall.cf) 
yield$root
```

### EXTRA

##### plot the yield curves (both maturity and call) if the bond is traded between 900 to 1200 THB

here are important part of code, please complete the rest

```{r,opts.label="badcode",class.source="badcodeStyle"}
  ##-- This is incomplete code block. It is UNABLE TO EXECUTE 

  tradeVal.List <- seq(900,1200,10)
  yield.List    <- ???      ##-- complete this line

  plot(tradeVal.List,yield.List)  
```

```{r,opts.label="hidden",class.source="hiddenStyle"}
  ##-- This is correct one

  tradeVal.List <- seq(900,1200,10)
  nList         <- length(tradeVal.List)
  yield.List    <- rep(0,length(tradeVal.List))
  for(i in 1:nList){
    bondCall.cf[length(bondCall.cf)] <- tradeVal.List[i]  
    yield.List[i]    <- uniroot(findNPV,interval=c(0,1),cashflow=bondCall.cf) 
  }
  
  plot(tradeVal.List,yield.List)  
```

------------------------------------------------------------------------

## Question 5 \## `r if (params$intab) "{.tabset}"`

##### Consider a 12-Years 7.980000% of annual coupon bond of BANPU plc issued on 01-Apr-11 with a par value of 1000 THB. Answer the following questions

------------------------------------------------------------------------

### 01: Yeild

##### a) If such bond trades at 1020 THB on 01-Apr-11, should we buy it?

```{r}
  banpu.cf <- c(-1020,rep(1000*0.0798,11),1000*1.0798)
  yeild <- uniroot(findNPV,interval=c(0,1),cashflow=banpu.cf) 
  yeild.banpu <- yeild$root
```

------------------------------------------------------------------------

### 02: Spread

##### b) What is the \textbf{spread} of this bond comparing to

[Thai Bond](http://ie.eng.chula.ac.th/~oran/classes/CompMeth_resources/scripts/thYeild-20130401.csv)

```{r,opts.label="badcode",class.source="badcodeStyle"}
  ##-- This code block is NOT EXECUTED

  ##-- download file first
  thaiYield <- read.csv(file="thYeild-20130401.csv",header=T)

  ##-- convert after using RStudio
  thaiYield  <- `thYeild-20130401` 
  mapThYield <- c(0,thaiYield$yield[4:15])/100
  npvSpread   <- function(cashflow,govYieldCurve,spread){
    ## cashflow <- banpu.cf
    ## govYieldCurve <-  mapThYield
    ## spread <- 0.0
    periods <- 1:length(cashflow) - 1
    npv     <- sum(cashflow*1/(1+govYieldCurve+spread)^periods )
  }
  root <- uniroot(npvSpread,c(0,0.2),cashflow=banpu.cf,govYieldCurve=mapThYield)
  spread <- root$root
```

#### EXTRA

##### repeat the Question 4.2 with current [Current Thai Bond](https://www.thaibma.or.th/EN/Market/YieldCurve/Government.aspx), assuming the other factors remains unchanged.

##### **hint:** here are some steps

-   import an excel file
-   remove unwanted columns
-   scale interesting columns, i.e., TTM (Time-To-Maturity) and Avg. Bidding Yield
-   interpolate/extrapolate the right TTM

```{r,opts.label="hidden",class.source="hiddenStyle"}
  ##-- This is correct one

  tradeVal.List <- seq(900,1200,10)
  nList         <- length(tradeVal.List)
  yeild.List    <- rep(0,length(tradeVal.List))
  for(i in 1:nList){
    bondCall.cf[length(bondCall.cf)] <- tradeVal.List[i]  
    yeild.List[i]    <- uniroot(findNPV,interval=c(0,1),cashflow=bondCall.cf) 
  }
  
  plot(tradeVal.List,yeild.List)  
```

------------------------------------------------------------------------

## Question 6 \## `r if (params$intab) "{.tabset}"`

##### Determine the ﬁrst- and second-order derivative at a specific point

------------------------------------------------------------------------

### 01: $f(x_1)$

##### a) $f(x_1) = sin(x_1)\,e^{-x_1^2}$, at $x = 1$

-   setup the function as expression

```{r}
  fExpr <- expression( sin(x1)*exp(-x1^2) )

  fExpr
  
  ##-- diff with respect to 'x1'
  D(fExpr,"x1")
```

-   evaluate diff of $f(x_1)$ single point at $x_1 = 1$, i.e.,

$$
  \nabla \mathbf{f}(\mathbf{x}) \equiv \frac{\partial}{\partial x_1} f(x_1 = 1)
$$

```{r}

  ##-- diff and evaluate 
  eval(D(fExpr,"x1"),list(x1=1))
```

-   use `function()` to evaluate many points

```{r}
  ##-- create dummy function
  fFirst <- function(x1){} 
  body(fFirst) <- D(fExpr,"x1")
  fFirst(1)
  
  fFirst(1:6)
  
```

-   evaluate diff of diff of $f(x_1)$ single point at $x_1 = 1$, i.e.,

$$
  \nabla^2 \mathbf{f}(\mathbf{x}) \equiv \frac{\partial^2}{\partial x_1^2} f(x_1 = 1)
$$

```{r}
  ##-- diff and evaluate 
  eval(D(D(fExpr,"x1"),"x1"),list(x1=1))
  
  ##-- create dummy function
  fSecond <- function(x1){} 
  body(fSecond) <- D(D(fExpr,"x1"),"x1")
  fSecond(1)
```

-   A common trick to evaluate diff of function at any point

```{r}
  ##-- common trick using `deriv3()` and create dummy function
  fAll <- function(x1){} 
  
  ##-- apply `derive3()`
  body(fAll) <- deriv3(fExpr,"x1")
  
  ##-- retrive value
  fAll(1)
  grad <-attr(fAll(1),"gradient")
  hess <-attr(fAll(1),"hessian")
```

### 02: $g(x_1,x_2)$

##### b) $g(x_1,x_2) = x_1\,e^{-x_2} + x_2 +1$, $(x_1,x_2) =(1,0)$

```{r,opts.label="badcode",class.source="badcodeStyle"}
  ##-- This code block is NOT EXECUTED
  gExpr <- expression( x1* exp(-1*x2) + x2 +1 )

  ##-- create dummy function
  gAll <- function(x1,x2){} 

  body(gAll) <- deriv3(gExpr,c("x1","x2"))
  gAll(1,0)
```

##### retrive $\nabla^2 \mathbf{g}(\mathbf{x})$ as matrix

```{r,opts.label="badcode",class.source="badcodeStyle"}
  ##-- This code is NOT execute
  matrix(attr(gAll(1,0),"hessian"),nrow=2)
```

------------------------------------------------------------------------

### EXTRA I

##### evaluate the first- and second derivative of $h(x_1,x_2,x_3,x_4)$

$$ 
h(x_1,x_2,x_3,x_4) = 100 (x_1-x_2)^2 + 10 (x_2-x_3)^4 + (x_3-x_4)^6 
$$

at $(x_1=1,x_2=2,x_3=-1,x_4=0)$

```{r,opts.label="badcode",class.source="badcodeStyle"}
  
  hExpr <-  expression( 100*(x1-x2)^2 + 10*(x2-x3)^4 + (x3-x4)^6 )
  hFunc <-  function(x1,x2,x3,x4){}
  body(hFunc) <- deriv3(hExpr,c("x1","x2","x3","x4"))

  hFunc(x1=1,x2=2,x3=-1,x4=0)[0]
  
  as.vector(attr(hFunc(x1=1,x2=2,x3=-1,x4=0),"gradient"))
  
  matrix(attr(hFunc(x1=1,x2=2,x3=-1,x4=0),"hessian"),nrow=4)
  
```

### EXTRA II

##### graph values of function, first-derivative, and second-derivative across their ranges

## Question 7 \##`r if (params$intab) "{.tabset}"`

##### If a company charges a price $p$ USD for a product, then it can sell $3000e^{-p}$ units of the product.

------------------------------------------------------------------------

### 01: Plot

##### a) plot revenue function

```{r}
### SOLUTION TO QUESTION 6A ###
 revFn    <- function(p){ 3000*p*exp(-1*p) }
 pRange   <- seq(0,5,0.1)
 revRange <- revFn(pRange)
 plot(pRange,revRange,cex=0.5,type="b",pch=16)
```

##### **Note**:

-   try `curve()` to smooth the function

### 02: Newton

##### b) use Newton method to find the optimal $p$, denoted by $p^*$, using $p=0$ as initial point for 10 iterations

-   setup expression

```{r}
### SOLUTION TO QUESTION 6B ### 
 revExp   <- expression(3000*p*exp(-1*p))
 rev1Deriv<- D(revExp,"p")
 rev2Deriv<- D(D(revExp,"p"),"p")
```

-   setup functions

```{r}
 rev1Val  <- function(p){}
 body(rev1Val) <- rev1Deriv 
 rev2Val  <- function(p){}
 body(rev2Val) <- rev2Deriv 
 revStep  <- function(p){-1*rev1Val(p)/rev2Val(p)  }
```

-   inital and execute loop

```{r}
 curP    <- 0
 for(i in 1:10){
   curObj  <- revFn(curP)
   curStep <- revStep(curP)
   nxtP    <- curP + curStep
   print(c(i,curP,curStep,curObj))
   curP    <- nxtP
 }
 
```

### 03: optimize()

##### c) use R function `optimize()` to solve for $p^*$

```{r}
### SOLUTION TO QUESTION 6C ### 
  pStar <- optimize(revFn,c(0,5))
  pStar$minimum
  pStar$objective
```

#### EXTRA

##### implement [Stochastic Search](https://www.cs.nott.ac.uk/~pszgxk/courses/g5aiai/003blindsearches/blind_searches.htm) and to find $p^*$

## Question 8 \## `r if (params$intab) "{.tabset}"`

##### Let $f(x) = 3\, x\, e^{-x^2} - 3 \frac{1}{x} \ln(x), x \geq\ 0$

### 01: Plot

##### a) plot $f(x)$ versus $x$ over the interval $[0,20]$

```{r}
### SOLUTION TO QUESTION 7A ###   
  fExpr <- expression(3*x*exp(-x^2)- 3*log(x)/x )  
  fFunc <- function(x){}
  body(fFunc) <- fExpr
  fDomain <- c(seq(0,1,0.01),seq(1,20,0.5) )
  fRange  <- sapply(fDomain,fFunc)
  plot(fDomain,fRange,cex=0.5,pch=16)
```

### 02: Newton

##### b) use the **Newton's Method** 10 iterations to locate $x_*$

-   setup expression

```{r}
### SOLUTION TO QUESTION 7A ###   
  fAll <- function(x){}
  body(fAll) <- deriv3(fExpr,"x")
```

-   initial and run loop

```{r}
  curX <- 2
  for(i in 1:10){
    ##i <- 1
    curX.all <- fAll(curX)
    curVal   <- curX.all[1]
    curGrad  <- attr(curX.all,"gradient")[1]
    curHess  <- attr(curX.all,"hessian")[1]    
    curStep  <- ifelse(curHess>0,-1*curGrad/curHess,runif(1))
    nextX    <- curX + curStep
    print(c(i,curX,curStep,curVal))
    ## browser() ##-- debug command
    curX    <- nextX
  }

```

-   compare result

```{r}
  optimize(fFunc,c(0,20),maximum=F)
```

# <!----------- BEGIN: STANDARD OKC FOOTER  -------------------------------->

------------------------------------------------------------------------

Copyright 2019   Oran Kittithreerapronchai.   All Rights Reserved.   **Last modified:** `r format(Sys.time(), '%Y-%M-%d, ')` <!----------- END: STANDARD OKC FOOTER  ---------------------------------->

<!-------- BEGIN: STANDARD OKC PRINT OUT with SYNTEX  -------------------->

```{r, echo=F,eval=F,include=F}

##--- CODE PRINING

curPath <- dirname(rstudioapi::getSourceEditorContext()$path)
fileName<- basename(rstudioapi::getSourceEditorContext()$path)
outName <- paste( c(fileName,".html"),collapse="")

setwd(curPath)

rmarkdown::render(fileName,params=list(intab=FALSE),output_file = outName)

# library(knitr)
# stitch(
#   script = fileName,
#   system.file("misc", "knitr-template.Rhtml", package = "knitr"),
#   output = outName
# )

##---
```

<!--------- ENDED: STANDARD OKC PRINT OUT with SYNTEX  ------------------->
