rm(list=ls()) # clears out memory

library(moments)

usd.data="./QUANDL_USDRUB.csv"
oil.data="./QUANDL_OPECOIL.csv"

url_usd="https://www.quandl.com/api/v1/datasets/CURRFX/USDRUB.csv?trim_start=2014-08-01"
download.file(url_usd, usd.data, "wget")

url_oil="https://www.quandl.com/api/v1/datasets/OPEC/ORB.csv?trim_start=2014-08-01"
download.file(url_oil, oil.data, "wget")

usd.data <- read.csv(usd.data, header = TRUE, sep = ",")
oil.data <- read.csv(oil.data, header = TRUE, sep = ",")

merged.data <- merge(usd.data,oil.data, by="Date")

usd_vs_oil = data.frame(DATE=merged.data$Date, 
                        USD=merged.data$Rate,
                        OIL=merged.data$Value)

remove(merged.data)

outliers <- c(97,112,127)

usd_vs_oil_clean = usd_vs_oil[-outliers,]


Ntot=dim(usd_vs_oil)[1]

LM1=lm(USD ~ I(OIL)+I(OIL^2),data=usd_vs_oil_clean)

print(summary(LM1)) #summary of results                                                                                                                       

## 
## Call:
## lm(formula = USD ~ I(OIL) + I(OIL^2), data = usd_vs_oil_clean)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1744 -1.2589 -0.1749  0.5070  9.0850 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.042e+02  3.701e+00  28.149  < 2e-16 ***
## I(OIL)      -1.024e+00  1.073e-01  -9.543  < 2e-16 ***
## I(OIL^2)     3.426e-03  7.312e-04   4.686 7.03e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Residual standard error: 2.637 on 128 degrees of freedom
## Multiple R-squared: 0.9427,  Adjusted R-squared: 0.9418 
## F-statistic:  1052 on 2 and 128 DF,  p-value: < 2.2e-16

## create plots

hist(LM1$residuals,nclass=24, freq=FALSE, col="lightgreen",
     main="Histogram of Residuals", xlab="Residuals")
lines(density(LM1$residuals),lwd=5, col="darkblue")

par(mfrow=c(1,1))
plot(LM1, which=c(4))

par(mfrow=c(2,2))    #drawing in 2 by 2 format                                                                                                                
plot(LM1)

par(mfrow=c(1,1))

tfirst = usd_vs_oil$DATE[1]
tlast = usd_vs_oil$DATE[Ntot]
ylast = usd_vs_oil$USD[Ntot]
xlast = usd_vs_oil$OIL[Ntot]

xmin = min(usd_vs_oil$OIL)
xmax = max(usd_vs_oil$OIL)

xx <- seq(xmin,xmax,len=200)
yy <- LM1$coef %*% rbind(1,xx,xx^2)


plot(usd_vs_oil$OIL,usd_vs_oil$USD,
     type="p",col="blue",lwd=4, xlab="OIL", ylab="USD",
     main=sprintf("Historical prices from  %s  to  %s",tfirst,tlast))

grid();
lines(xx,yy,lwd=4,col="magenta")
points(xlast, ylast, type = "p", col="red", lwd=5, pch = 21)
text(xlast, ylast, tlast, adj=-0.2, col="darkred")

cat("mean=",mean(LM1$resid),"\n")

## mean= 5.198874e-17

cat("st.dev.=",sd(LM1$resid),"\n")

## st.dev.= 2.616526

cat("skewness=",skewness(LM1$resid),"\n")

## skewness= 1.54103

cat("kurtosis=",kurtosis(LM1$resid),"\n")

## kurtosis= 6.536955