# PREMIUM 

premium<-10

# TIME HORIZON 

tm<-60

p<-matrix(c(rep(premium,tm/2),rep(0,tm/2)),tm,1)

# NUMBER OF SIMULATIONS

nreplic<-10000

# GENERATE FROM THE LOGNORMAL DISTRUBUTION THE INTEREST RATES

mean<-0  #SET THE MEAN EQUAL TO ZERO FOR Q-MEASURE

std<-0.155

m<-log(mean+1)

meanvec=m*matrix(1,tm,1)
stdvec=std*matrix(1,tm,1)

yr1=matrix(0,nreplic,tm)

stockbond <- function(p, rates, nrep, y, porc) {

s<-matrix(0,tm,1)
b<-matrix(0,tm,1)
r<-matrix(0,tm,1)
t<-matrix(0,tm,1)
f<-matrix(0,tm,1)
e<-matrix(0,tm,1)
y<-matrix(0,tm,1)

qorc=1-porc;
 s[1,1]=p[1,1]*porc[1,1]
 b[1,1]=p[1,1]*qorc[1,1]

 r[1,1]=0
 t[1,1]=s[1,1]+b[1,1]
 f[1,1]=s[1,1]*(1+rates[1,1])
 e[1,1]=b[1,1]+f[1,1]

    for (i in 2:tm){
        if (i<=tm/2) {
           r[i,1]<-0
           }
        if (i>tm/2)  {
           r[i,1]<-premium
           }
        t[i,1]<-e[i-1,1]+p[i,1]-r[i,1]
        if (t[i,1]>0){
           s[i,1]<-t[i,1]*porc[i,1]
           b[i,1]<-t[i,1]*qorc[i,1]
           f[i,1]<-s[i,1]*(1+rates[i,1])
           e[i,1]<-b[i,1]+f[i,1]
           }
        if (t[i]<=0) {
           e[i,1]<-t[i,1]
           }
     }
 return(e)
 }

alfa<-0.2
sss1<-0.00
sss2<-0.01

timepension <- function(p, rates, nrep, y, alfa,sss) {

s<-matrix(0,tm,1)
b<-matrix(0,tm,1)
r<-matrix(0,tm,1)
t<-matrix(0,tm,1)
f<-matrix(0,tm,1)
e<-matrix(0,tm,1)
y<-matrix(0,tm,1)

c<-matrix(0,tm,1)
d<-matrix(0,tm,1)
u<-matrix(0,tm,1)

c[1,1]<-p[1,1]
d[1,1]<-p[1,1] 
u[1,1]<-0 


 s[1,1]<-p[1,1]
 b[1,1]<-0

 r[1,1]<-0

 for (i in 2:tm){

  if (i<=30) {
      r[i,1]<-0
      }  
  if (i>30) {
      r[i,1]<-10
      }  
  if (c[i-1,1]>0) {
     f[i-1,1]<-c[i-1,1]*(1+rates[i-1,1])
     c[i,1]<-f[i-1,1]+p[i,1]-r[i,1]
     d[i,1]<-d[i-1,1]*(1+sss)+p[i,1]-r[i,1]+alfa*(c[i-1,1]*(1+rates[i,1]) -(1+sss)*d[i-1,1])
     u[i,1]<-(1-alfa)*(c[i-1,1]*(1+rates[i,1])-d[i-1,1]*(1+sss))
     }
  if (c[i-1,1]<=0) {
     c[i,1]<-c[i-1,1]+p[i,1]-r[i,1]
     if (d[i-1,1]>0) { 
        d[i,1]<-d[i-1,1]*(1+sss)+p[i,1]-r[i,1]+alfa*(c[i-1,1]*(1) -(1+sss)*d[i-1,1])
        u[i,1]<-(1-alfa)*(c[i-1,1]*(1)-d[i-1,1]*(1+sss))
        }
     if (d[i-1,1]<=0) { 
        d[i,1]<-d[i-1,1]*(1)+p[i,1]-r[i,1]+alfa*(c[i-1,1]*(1) -(1)*d[i-1,1])
        u[i,1]<-(1-alfa)*(c[i-1,1]*(1)-d[i-1,1]*(1))      
        }  
     }
  }
 return(d)
}


# SIMULATION STARTS


set.seed(1)

for (nrep in 1:nreplic){

  z<-matrix(0,tm,1)
  z[,1]<-rnorm(tm,0,1)
  lns<-(z*stdvec)+meanvec
  rates<-exp(lns)-1


# IN THIS PART OF THE PROGRAM ONE CAN CHOOSE THE PENSION PRODUCT
# YOU CAN ALSO SPECIFY ANY PENSION PRODUCT IN porc<- matrix(...)

  # Unit Linked
   porc<-matrix(0.5,tm,1)  

  # Stocks 50
  # porc<-matrix(c(rep(0.5,tm/2),rep(0,tm/2)),tm,1) 

  # Market Pension
  # porc<-matrix(c(rep(0.6,13),rep(0.5,8),rep(0.4,6),rep(0.3,4),rep(0.2,5),rep(0.1,5),rep(0,19)),tm,1)    

  # Danica
  # porcaux<-c(0.750, 0.729, 0.674, 0.626, 0.586, 0.551, 0.520, 0.493, 0.469,0.447, 
  # 0.428, 0.410, 0.394, 0.380, 0.366, 0.354, 0.343, 0.333, 0.323, 0.314, 0.306,
  # 0.298, 0.290, 0.284, 0.277, 0.271, 0.265,0.260, 0.255,0.250)
  # edadindex<-seq(from=60,to=79,by=1)
  # porcaux2<-0.25-0.25*(edadindex-59)/20
  # porcaux3<-rep(0,10)
  # porc2<-c(porcaux,porcaux2,porcaux3)
  # porc<-matrix(porc2,tm,1) 

  # Nordea
  # edadindex=seq(from=30,to=89,by=1)
  # porcaux=0.44-(0.37/60)*(edadindex-30) 
  # porc<-matrix(porcaux,tm,1)

  # PFA
  # porc<-matrix(c(rep(0.77,25),rep(0.52,10),rep(0.27,25)),tm,1) 

   yr1[nrep,]<-stockbond(p,rates,nrep,y,porc)

  # TIMEPENSION
  # newrates<-0.7*rates

  # yr1[nrep,]<-timepension(p, newrates, nrep, y, alfa, sss1)
    }


# Table 5 Fair price of a guarantee

mean<-0

m0<-log(mean+1)

meanvec0=m0*matrix(1,tm,1)

# SIMULATION STARTS FOR PRICE GUARANTEE NON PATH DEPENDENT PRODUCTS

yr1pg=matrix(0,nreplic,tm)

set.seed(1)

for (nrep in 1:nreplic){

  z<-matrix(0,tm,1)
  z[,1]<-rnorm(tm,0,1)
  lns0<-(z*stdvec)+meanvec0
  rates0<-exp(lns0)-1
  yr1pg[nrep,]<-stockbond(p,rates0,nrep,y,porc)
  }

 b<-yr1pg[,60][yr1pg[,60]<0]
 price_guar=-(mean(b)*(length(b)/10000))
 price_guar

# SIMULATION STARTS FOR PRICE GUARANTEE TIMEPENSION

mean<-0

m0<-log(mean+1)

meanvec0=m0*matrix(1,tm,1)

yr1pg=matrix(0,nreplic,tm)

set.seed(1)

for (nrep in 1:nreplic){

  z<-matrix(0,tm,1)
  z[,1]<-rnorm(tm,0,1)
  lns0<-(z*stdvec)+meanvec0
  rates0<-exp(lns0)-1
  newrates0<-0.7*rates0
  yr1pg[nrep,]<-timepension(p, newrates0, nrep, y, alfa, sss1)
  }

 b<-yr1pg[,60][yr1pg[,60]<0]
 price_guar=-(mean(b)*(length(b)/10000))
 price_guar