############################
###    VEZBE IZ EFM-A    ###
############################

# install.packages("FinCal")
library(FinCal)

####################
## PRVO POGLAVLJE ##
####################

## Funkcije u R-u (paket FinCal): 
      # npv(r,cf)
      # irr(cf)
      # r.continuos(r,m)
      # r.norminal(rc, m)

## 1.1
# unutrasnja stopa dobiti
irr(cf = c(-100000,0,52000,0,54000)) 
# sadasnja vrednost za r=0.1
npv(r = 0.1,cf = c(-100000,0,52000,0,54000))

## 1.2. 
npv(r = 0.1,cf = rep(500000,20))

## 1.3.
# za sest mecesi
(pv1 <- npv(r = 0.12, rep(1000,6)))
(pv2 <- npv(r = 0.12, c(1000,rep(900,6))))
# za godinu dana
(pv1 <- npv(r = 0.12, rep(1000,12)))
(pv2 <- npv(r = 0.12, c(1000,rep(900,12))))

## 1.4. 
500000*(1.07/1.03)^10

## 1.7.
2000000*(1+0.015)/(200*(1+0.01))

## 1.8.
70000/(sum((1+0.01/12)^(1:(12*16))))

## 1.9.
1000*(1+0.02/2)^2*(1+0.03/2)^2


#####################
## DRUGO POGLAVLJE ##
#####################

## 2.2.
# deo a)
N <- 10
t <- c(0.5,1,1.5,1.75,2,2.5,3.1,4,5,6) 
B <- c(9.75,9.5,9.1,8.5,8.3,8,7.58,7.2,7,6.8)
(R <- (N/B)^(1/t)-1)
# deo b)
model1 <- lm(R ~ t+I(t^2)+I(t^3))
summary(model1)
model2 <- lm(R ~ t+I(t^2))
summary(model2)
# p vrednost je manja za model2
anova(model1,model2) # velika p vrednost, pa prihvatamo H0 i biramo jednostavniji model 
plot(t,R,ylim = c(0,0.1))
lines(t,0.01513+0.064846*t-0.016273*t^2+0.001143*t^3, col='deeppink')
lines(t,0.033631+0.036296*t-0.005335*t^2, col='darkturquoise')
# deo c)
newdata1 <- data.frame(t=c(8/12,16/12,24/12,32/12))
x <- predict(model2,newdata1)+1
(PV1 <-100 - 20/(x[1]^(2/3))+50/(x[2]^(4/3))+10/(x[3]^2)+20/(x[4]^(8/3)))
(PV2 <- -20 + 100/(x[1]^(2/3))+10/(x[2]^(4/3))+30/(x[3]^2))

## 2.3. 
tok <- c(-10100,1000,1000,11000)
irr(tok) # unutrasnja stopa dobiti

## 2.4. 
# deo a)
a <- 1/1.0495
(fer_cena <- 5*sum(a^(c(1,2,3,4,5)))+100*(a^5))
# ili
npv(r = 0.0495, cf = c(0,5,5,5,5,105))
# deo b)
x <- seq(0,0.3,0.05)
y <- c()
for(i in 1:length(x))
y[i] <- 5*sum( 1/((x[i]+1)^(1:5))) + 100*(1/((x[i]+1)^5))
plot(x,y, type = "l")
points(0.0495,100.2168)

## 2.6.
t <- seq(0.5,5,0.5)
r <- 0.05-0.03*(1-exp(-t))/t 
250*sum(1/((1+r)^t))+10000/((1+r[10])^5)

## 2.7.
# install.packages("YieldCurve")
library(YieldCurve)
# podaci iz zadatka 2.2
t <- c(0.5,1,1.5,1.75,2,2.5,3.1,4,5,6) 
B <- c(9.75,9.5,9.1,8.5,8.3,8,7.58,7.2,7,6.8)
(R <- (N/B)^(1/t)-1)
ModelNS <- Nelson.Siegel(R,t)
ModelNS
# install.packages("NMOF")
library(NMOF)
ModelOcena <- NS(ModelNS[1,],t)
ModelOcena 
# grafik
plot(t,R,ylim = c(0,0.1))
lines(t,0.01513+0.064846*t-0.016273*t^2+0.001143*t^3, col='deeppink')
lines(t,0.033631+0.036296*t-0.005335*t^2, col='darkturquoise')
lines(t,ModelOcena, col='green')

## 2.17.
delta_t<-1/2
R <- c(3,3.5,3.75,4,5.1,5.25)/100
t <- seq(0,2.5,0.5)
(f<-((1+R[2*(t+delta_t)])^(t+delta_t)/(1+R[2*delta_t])^delta_t)^(1/t)-1)
(P <- sum(4/(1+R)^(t+1/2))+100/((1+R[6])^(t[6]+1/2)))


#####################
## TRECE POGLAVLJE ##
#####################

## 3.1.
S <- matrix(c(1,0,0.5,0,2,0.5,0.5,0.5,1),3,3)
S_inverz <- solve(S)
(omega <- rowSums(S_inverz)/sum(S_inverz))
(r_min <- sum(omega*c(0.1,0.15,0.1)))
(S_min <- t(omega)%*%S%*%omega)

## 3.2.
S <- matrix(c(0.07,-0.005,-0.005,0.013),2,2)
S_inverz <- solve(S)
R <- c(0.175,0.055)
# portfolio minimalne disperzije
(omega <- rowSums(S_inverz)/sum(S_inverz))
100000*omega
(r_min <- sum(omega*R))
(S_min <- t(omega)%*%S%*%omega)
# VaR
(VaR <- 100000*(2.33*sqrt(S_min)-r_min)) 
# grafik
omega <- seq(from=-0.5,to=1.5,by=0.05)
sigma <- sqrt(S[1,1]*omega^2+S[2,2]*(1-omega)^2+2*S[1,2]*omega*(1-omega))
r <- omega*R[1]+(1-omega)*R[2]
plot(sigma,r,ylim=c(0,0.3),xlim=c(0,0.28),type="l",xlab=expression(sigma),lty=2)
legend('topright',"efkasni portfolio",col="deeppink3",lty=1,cex=0.6)
sigm_E <- sigma[which.min(sigma):length(sigma)]
r_E <- r[which.min(sigma):length(sigma)]
lines(sigm_E,r_E,ylim=c(0,0.3),type="l",
      xlab=expression(sigma),col="deeppink3",lwd=2)
# tangentni portflio
rf <- 0.005
(omega_t <- S_inverz%*%(R-rf)/sum(S_inverz%*%(R-rf)))
(r_t <- sum(omega_t*R))
(S_t <- t(omega_t)%*%S%*%omega_t)
# grafik - Efikasni portfolio koji se sastoji od akcija i obveznice
abline(rf,(r_t-rf)/sqrt(S_t))
points(sqrt(S[2,2]),R[2],col='orange',lwd=2)
points(sqrt(S[1,1]),R[1],col='orange',lwd=2)
points(sqrt(S_min),r_min,col='green',lwd=2)
points(0,rf,col='blue',lwd=2)
points(sqrt(S_t),r_t,col='blue',lwd=2)
points(0.059,0.055,col='red',lwd=2)

## 3.3.
R <-c(0.0427,0.0016,0.0285)
S <- matrix(c(0.01,0.0018,0.0011,0.0018,0.0109,0.0026,0.0011,0.0026,0.0199),3,3)
rf <- 0.005
S_inverz <- solve(S)
(omega_t <- S_inverz%*%(R-rf)/sum(S_inverz%*%(R-rf)))
(S_t <- t(omega_t)%*%S%*%omega_t)
(S_3 <- t(c(1/3,1/3,1/3))%*%S%*%c(1/3,1/3,1/3))
w_a <- seq(from=-0.5,to=1.5,by=0.05)
w_b <- seq(from=-0.5,to=1.5,by=0.05)
r <- matrix(0, length(w_a), length(w_b))
sigma <- matrix(0, length(w_a), length(w_b))
for(i in 1:length(w_a)){
  for(j in 1:(length(w_b))){
    r[i,j] <- R[1]*w_a[i]+R[2]*w_b[j]+R[3]*(1-w_a[i]-w_b[j])
    sigma[i,j] <- sqrt(S[1,1]*w_a[i]^2+S[2,2]*w_b[j]^2+S[3,3]*(1-w_a[i]-w_b[j])^2
                       + 2*S[1,2]*w_a[i]*w_b[j]+2*S[1,3]*w_a[i]*(1-w_a[i]-w_b[j])+ 2* S[2,3]* w_b[j]*(1-w_a[i]-w_b[j]))
  }
}
plot(sigma,r, xlim = c(0.07,0.2), ylim = c(-0.01,0.06))

## 3.8.
# portfolio minimalne disperzije 
S <- matrix(c(0.0064,-0.0024,-0.0024,0.01),2,2)
S_inverz <- solve(S)
R <- c(0.05,0.12)
(omega <- rowSums(S_inverz)/sum(S_inverz))
(r_min <- sum(omega*R))
(S_min <- t(omega)%*%S%*%omega)
# tangentni portflio
rf <- 0.01
(omega_t <- S_inverz%*%(R-rf)/sum(S_inverz%*%(R-rf)))
(r_t <- sum(omega_t*R))
(S_t <- t(omega_t)%*%S%*%omega_t)
(Sarpov_kolicnik <- (r_t-rf)/sqrt(S_t))


####################
## PETO POGLAVLJE ##
####################

## 5.1-5.5 Binomni model
BinomniModel <- function(TypeFlag=c("CE","PE","CA","PA"),S,K,t,r,d,u,n){
  
  if(TypeFlag=="CE" || TypeFlag=="CA") z <- 1
  if(TypeFlag=="PE" || TypeFlag=="PA") z <- -1

  dt <- t/n
  p <- (exp(r*dt)-d)/(u-d)
  Df <- exp(-r*dt)
  OptionValue <- z*(S*u^(0:n)*d^(n:0)-K)
  OptionValue <- (abs(OptionValue)+OptionValue)/2
  
  if(TypeFlag=="CE" || TypeFlag=="PE"){
    for (j in (n-1):0) 
      for (i in 0:j)
        OptionValue[i+1] <- (p*OptionValue[i+2]+(1-p)*OptionValue[i+1])*Df
  }
  
  if(TypeFlag=="CA" || TypeFlag=="PA"){
    for (j in (n-1):0) 
      for (i in 0:j)
        OptionValue[i+1] <- max((z*(S*u^i*d^(abs(i-j))-K)),(p*OptionValue[i+2]+(1-p)*OptionValue[i+1])*Df) 
  }
  return(OptionValue)
}
## 5.1.
BinomniModel(TypeFlag="CE",S=50,K=50,t=0.5,r=0.01,u=55/50,d=45/50,n=1)
BinomniModel(TypeFlag="PE",S=50,K=50,t=0.5,r=0.01,u=55/50,d=45/50,n=1)
## 5.2.
BinomniModel(TypeFlag="CE",S=45,K=46,t=1/3,r=0.1,d=0.92,u=1.1,n=2)
## 5.3.
BinomniModel(TypeFlag="PE",S=50,K=46,t=1/3,r=0.05,u=1.09,d=0.9,n=2)
BinomniModel(TypeFlag="PA",S=50,K=46,t=1/3,r=0.05,u=1.09,d=0.9,n=2)
## 5.4.
BinomniModel(TypeFlag="PA",S=45,K=46,t=1,r=0.1,u=1.1,d=0.9,n=3)
## 5.5.
BinomniModel(TypeFlag="CE",S=50,K=50,t=1,r=0.05,1.05,0.95,n=4)
BinomniModel(TypeFlag="PA",S=50,K=50,t=1,r=0.05,u=1.05,d=0.95,n=4)

## 5.9.
(PV <- exp(-0.1/12)+exp(-0.1/6)+exp(-0.1/3))
50-PV
(d1 <- (log(47.058/50)+(0.1+0.09/2)/2)/(0.3*sqrt(1/2)))
(d2 <- d1-(0.3*sqrt(1/2)))
(p <- 50*exp(-0.1/2)*pnorm(0.156)-50*pnorm(-0.056))

## 5.12.
BS <-   function(S, K, t, r, sig, type="C"){
  d1 <- (log(S/K) + (r + sig^2/2)*t) / (sig*sqrt(t))
  d2 <- d1 - sig*sqrt(t)
  if(type=="C"){
    value <- S*pnorm(d1) - K*exp(-r*t)*pnorm(d2)
  }
  if(type=="P"){
    value <- K*exp(-r*t)*pnorm(-d2) - S*pnorm(-d1)
  }
  return(value)
}
implied.vol <- function(S, K, t, r, market, type){
  sig <- 0.20
  sig.up <- 1
  sig.down <- 0.001
  err <- BS(S, K, t, r, sig, type) - market 
  while(abs(err) > 0.00001){
    if(err < 0){
      sig.down <- sig
      sig <- (sig.up + sig)/2
    }else{
      sig.up <- sig
      sig <- (sig.down + sig)/2
    }
    err <- BS(S, K, t, r, sig, type) - market
  }
  return(sig)
}
implied.vol(100,110,1/2,0.05,10,type="C")