PBSR assignment rmd.Rmd

---
title: "PBSR Assignment 2"
author: "Sayantika Sengupta, Aniket Saha, Anjan"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Problem 3:

```{r}
library(scales)
# library(AICcmodavg)
```

```{r}
attach(faithful)
hist(faithful$waiting,xlab = 'waiting',probability = T,col='pink',main='')

```

```{r}
data_q3 <- faithful 
x <- sort(data_q3$waiting)
```

comparing 3 models

```{r}
# model 1
p <- length(x[x<65])/length(x)
as <- mean(x[x<65])
ass <- var(x[x<65])
s <- ass/as
a <- as/s
mu <- mean(x[x>=65])
sigma <- sd(x[x>=65])
theta_inital <- c(p, a, s, mu, sigma)
neg_log_likelihood <- function(theta, data){
    n = length(data)
    
    p = theta[1]
    a = theta[2]
    s = theta[3]
    mu = theta[4]
    sigma = theta[5]
    
    l = 0
    for (i in 1:n) {
        l = l + log(p*dgamma(data[i], shape = a, scale = s) + (1-p)*dnorm(data[i], mean = mu, sd = sigma))
    }
    return(-l)
}
fit = optim(theta_inital,
      neg_log_likelihood,
      data = x,
      control = list(maxit = 1500),
      lower = c(0, 0, 0, -Inf, 0),
      upper = c(1, Inf, Inf, Inf, Inf),
      method="L-BFGS-B")
theta_1 = fit$par
theta_1
p = theta_1[1]
a = theta_1[2]
s = theta_1[3]
mu = theta_1[4]
sigma = theta_1[5]
model_1 = p*dgamma(x, shape = a, scale = s) + (1-p)*dnorm(x, mean = mu, sd = sigma)
aic_1 <- 2*5 + neg_log_likelihood(theta_1, x)

```

```{r}
# model 2
p <- length(x[x<65])/length(x)
as_1 <- mean(x[x<65])
ass_1 <- var(x[x<65])
s_1 <- ass_1/as_1
a_1 <- as_1/s_1
as_2 <- mean(x[x>=65])
ass_2 <- var(x[x>=65])
s_2 <- ass_2/as_2
a_2 <- as_2/s_2
theta_inital <- c(p, a_1, s_1, a_2, s_2)
neg_log_likelihood <- function(theta, data){
    n <- length(data)
    
    p <- theta[1]
    a_1 <- theta[2]
    s_1 <- theta[3]
    a_2 <- theta[4]
    s_2 <- theta[5]
    
    l <- 0
    for (i in 1:n) {
        l = l + log(p*dgamma(data[i], shape = a_1, scale = s_1) + (1-p)*dgamma(data[i], shape = a_2, scale = s_2))
    }
    return(-l)
}
fit = optim(theta_inital,
      neg_log_likelihood,
      data = x,
      control = list(maxit = 1500),
      lower = c(0, 0, 0, 0, 0),
      upper = c(1, Inf, Inf, Inf, Inf),
      method="L-BFGS-B")
theta_2 <- fit$par
theta_2
p <- theta_2[1]
a_1 <- theta_2[2]
s_1 <- theta_2[3]
a_2 <- theta_2[4]
s_2 <- theta_2[5]
model_2 <- p*dgamma(x, shape = a_1, scale = s_1) + (1-p)*dgamma(x, shape = a_2, scale = s_2)
aic_2 <- 2*5 + neg_log_likelihood(theta_2, x)

```

```{r}
# model 3
p <- length(x[x<65])/length(x)
m_1 <- mean(x[x<65])
v_1 <- var(x[x<65])
sigma2_1 <- log((v_1/m_1^2) + 1)
mu_1 <- log(m_1) - sigma2_1/2
m_2 <- mean(x[x>=65])
v_2 <- var(x[x>=65])
sigma2_2 <- log((v_2/m_2^2) + 1)
mu_2 <- log(m_2) - sigma2_2/2
theta_inital <- c(p, mu_1, sqrt(sigma2_1), mu_2, sqrt(sigma2_2))
neg_log_likelihood <- function(theta, data) {
    n <- length(data)
    
    p <- theta[1]
    mu_1 <- theta[2]
    sigma_1 <- theta[3]
    mu_2 <- theta[4]
    sigma_2 <- theta[5]
    
    l <- 0
    for (i in 1:n) {
        l = l + log(p*dlnorm(data[i], meanlog = mu_1, sdlog = sigma_1) + (1-p)*dlnorm(data[i], meanlog = mu_2, sdlog = sigma_2))
    }
    
    return(-l)
}
fit = optim(theta_inital,
      neg_log_likelihood,
      data = x,
      control = list(maxit = 1500),
      lower = c(0, -Inf, 0, -Inf, 0),
      upper = c(1, Inf, Inf, Inf, Inf),
      method="L-BFGS-B")
theta_3 <- fit$par
theta_3
p <- theta_3[1]
mu_1 <- theta_3[2]
sigma_1 <- theta_3[3]
mu_2 <- theta_3[4]
sigma_2 <- theta_3[5]
model_3 <- p*dlnorm(x, meanlog = mu_1, sdlog = sigma_1) + (1-p)*dlnorm(x, meanlog = mu_2, sdlog = sigma_2)
aic_3 <- 2*5 + neg_log_likelihood(theta_3, x)
```

```{r}
hist(x, xlab = 'waiting', probability = T, col='pink', main='')
lines(x, model_1, col = "red")
lines(x, model_2, col = "blue")
lines(x, model_3, col = "green")
```

```{r}
results <- data.frame(
    models = c("Gamma + Normal", "Gamma + Gamma", "Lognormal + Lognormal"),
    AIC = c(aic_1, aic_2, aic_3)
)
results
```

Based on the AIC value of all the models, we choose the third model as it has the lowest AIC value.

The required probability $\mathbb{P}[60<\texttt{waiting}<70]$ is:
```{r}
p <- theta_3[1]
mu_1 <- theta_3[2]
sigma_1 <- theta_3[3]
mu_2 <- theta_3[4]
sigma_2 <- theta_3[5]
reqd_prob <- (p*plnorm(70, meanlog = mu_1, sdlog = sigma_1) + (1-p)*plnorm(70, meanlog = mu_2, sdlog = sigma_2)) - (p*plnorm(60, meanlog = mu_1, sdlog = sigma_1) + (1-p)*plnorm(60, meanlog = mu_2, sdlog = sigma_2))
reqd_prob
```


## Problem 5:

```{r,warning=FALSE,message=FALSE}
library(quantmod)
getSymbols('TCS.NS')#getting the tcs data
tail(TCS.NS)
plot(TCS.NS$TCS.NS.Adjusted)#visualizing the tcs data
```
```{r,warning=FALSE,message=FALSE}
getSymbols('^NSEI')#getting nsei dataset
tail(NSEI)
plot(NSEI$NSEI.Adjusted)#plotting the nsei data
```
```{r,warning=FALSE,message=FALSE}
TCS_rt = diff(log(TCS.NS$TCS.NS.Adjusted))
Nifty_rt = diff(log(NSEI$NSEI.Adjusted))
retrn = cbind.xts(TCS_rt,Nifty_rt)
retrn = na.omit(data.frame(retrn))
plot(retrn$NSEI.Adjusted,retrn$TCS.NS.Adjusted
     ,pch=20
     ,xlab='Market Return'
     ,ylab='TCS Return'
     ,xlim=c(-0.18,0.18)
     ,ylim=c(-0.18,0.18))
grid(col='grey',lty=1)
```
It is given that 
$$
r_{tcs} = \alpha + \beta r_{ni} + \epsilon
$$
### MME method:
```{r,warning=FALSE,message=FALSE}
rtcs = mean(retrn$TCS.NS.Adjusted)
rni = mean(retrn$NSEI.Adjusted)
row = cor(retrn$TCS.NS.Adjusted,retrn$NSEI.Adjusted)
sdtcs = sqrt(var(retrn$TCS.NS.Adjusted))
sdnifty = sqrt(var(retrn$NSEI.Adjusted))
n = nrow(retrn)
b1 = row*(sdtcs/sdnifty)
a1 = sdtcs - b1*sdnifty
rtcs_hat = a1 + b1 * retrn$NSEI.Adjusted
error = retrn$TCS.NS.Adjusted - rtcs_hat
sigma1 = sqrt(var(error))
Method_of_Moments <- c(a1,b1,sigma1)
Method_of_Moments
```

### Ordinary least square method:

```{r,warning=FALSE,message=FALSE}
lin_mod = summary(lm(TCS.NS.Adjusted~NSEI.Adjusted, data = retrn))
a2 = lin_mod$coefficients [1,1]
b2 = lin_mod$coefficients [2,1]
Method_of_Moments <- c(a1,b1,sigma1)
k = (n-2)/n
errornew = retrn$TCS.NS.Adjusted - (a2+ b2*retrn$NSEI.Adjusted)
sigma2 = (sqrt(var(errornew)))*k
OLS <- c(a2,b2, sigma2)
```

### Representing the estimates as a dataframe:

```{r}
Parameters <- c("alpha", "beta", "sigma")
table = data.frame(Parameters, Method_of_Moments,OLS)
table
```

### Estimated rise in price of tcs

```{r}
est1 = 3200 - ( 200*b1)
est2 = 3200 - ( 200*b2)
est = c(est1,est2)
est
```