##########
#Practice 7
#
##########

#HW1 from last week is prolonged
#I suggested a solution based on simulating from the 
#joint (normal) distribution of the parameter estimators
#(fbvevd in the evd package provides this) and then the
#for the grid points it can be checked if they belong to 
#the region or not. But it is still far from obvious

#HW2
library(POT)
sizes = c(40,100,200,300)
rep=100
res=array(0,dim=c(rep,length(sizes),3))
p = 0.95
set.seed=1234
for (j in 1:rep)
{
for (i in 1:length(sizes)) {
  ech = rgpd(sizes[i] ,shape=0.6,scale=1) 
  q = quantile(ech,p) #mle estimation 
  fit = fitgpd(ech, threshold=0, est="mle") # mle estimation of parameters 
  pli = confint(fit, prob=p, parm = "quant", level=0.95,range=c(2,55), prof=TRUE) #prof. lik. int
  res[j,i,1] = pli[2]-pli[1] 
  res[j,i,2] = pli[1]
  res[j,i,3] = pli[2]
	}
}
rownames(res)=as.character(sizes)
colnames(res)=c("value","b_inf","b_sup")
print(res[,1,])
theo=qgpd(0.95,shape=0.6,scale=1)
cover=c(0,0,0,0)
for (i in 1:length(sizes)) {
	for (j in 1:rep)
	{if (!is.na(res[j,i,2]) & !is.na(res[j,i,3])) 
		{if (theo>res[j,i,2] & theo<res[j,i,3]) cover[i]=cover[i]+1}
	}
}
print(theo)
cover
#79 97 95 95 #but it might be higher if done manually and the NA's could 
#be eliminated


#for the ML
res=array(0,dim=c(rep,length(sizes),3))
p = 0.95
set.seed=1234
for (j in 1:rep)
{
for (i in 1:length(sizes)) {
  ech = rgpd(sizes[i] ,shape=0.6,scale=1) 
  q = quantile(ech,p) #mle estimation 
  fit = fitgpd(ech, threshold=0, est="mle") # mle estimation of parameters 
  pli = confint(fit, prob=p, parm = "quant", level=0.95, prof=FALSE) #ML. lik. int
  res[j,i,1] = pli[2]-pli[1] 
  res[j,i,2] = pli[1]
  res[j,i,3] = pli[2]
	}
}
rownames(res)=as.character(sizes)
colnames(res)=c("value","b_inf","b_sup")
print(res[,1,])
theo=qgpd(0.95,shape=0.6,scale=1)
cover=c(0,0,0,0)
for (i in 1:length(sizes)) {
	for (j in 1:rep)
	{if (!is.na(res[j,i,2]) & !is.na(res[j,i,3])) 
		{if (theo>res[j,i,2] & theo<res[j,i,3]) cover[i]=cover[i]+1}
	}
}
print(theo)
cover #89 95 88 90 so the prof lik is slightly better in this range


#bootstrap
#####
#task 1
#####

library(boot)

#A simple regression example
y = read.table("http://zempleni.elte.hu//nasdaq_djia_2005_2010.txt")

betfun = function(data,b,formula){  
# b is the random indices for the bootstrap sample
	d = data[b,] 
	return(lm(d[,1]~d[,2], data = d)$coef[2])  
# that is for the beta coefficient
	}
# now you can bootstrap:
bootbet = boot(data=y, statistic=betfun, R=500) 
# R is how many bootstrap samples
names(bootbet)
plot(bootbet)
hist(bootbet$t) 
 lm1=lm(y[,1]~y[,2])
error=summary(lm1)$coefficients[2,2]*1.96
left=summary(lm1)$coefficients[2,1]-error
right=summary(lm1)$coefficients[2,1]+error

#draw the normality-based confidence interval
abline(v=left,col=2,lty=2)
abline(v=right,col=2,lty=2)
#what is our opinion?

####
#Task 2
####

################
#bootstrap bca konf. int.
################
library(boot)
plot (cd4)
# boot sample quantiles
ii=c(1:dim(cd4)[1])
n=5000
ered=rep(0,times=n)
for (i in 1:n)
{ss=sample(ii,replace=T)
ered[i]=corr(cd4[ss,])
}
conf1=c(quantile(ered,0.005),quantile(ered,0.995))

# confidence intervals for the medical data
#ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
cd4.boot <- boot(cd4, corr, R = 999, stype = "w", sim = "ordinary")
boot.ci(cd4.boot, conf = c(0.95,0.99, 0.995),type = "bca")
#normal approximation (not valid)...
error=2.58*(1-corr(cd4)^2)/sqrt(dim(cd4)[1])

##
#Task 3: limitations of iid bootstrap
##
+
##########
#Task 3/2: AR(1) models
##########

gen=500  # bootstrap repetitions
nnn=5*10^3
library(np)

set.seed(12345)
alfa=0.1
yt <- arima.sim(n=nnn, list(ar=c(alfa)))
par(mfrow=c(2,2))
plot(yt[1:100],type="l")
plot(ARMAacf(ar=c(alfa), lag=50), type="h")
plot(ARMAacf(ar=c(alfa), lag=50, pacf=TRUE), type="h")
b.star(yt,round=TRUE)  # optimal block size, Politis-White algorithm

alfa=0.9
yt <- arima.sim(n=nnn, list(ar=c(alfa)))
par(mfrow=c(2,2))
plot(yt[1:100],type="l")
plot(ARMAacf(ar=c(alfa), lag=50), type="h")
plot(ARMAacf(ar=c(alfa), lag=50, pacf=TRUE), type="h")

b.star(yt,round=TRUE)  # optimal block size

# i.i.d. bootstrap
nn=gen
mea=rep(0,times=gen)
for (i in 1:nn)  mea[i]=mean(sample(yt,replace=T))
becs1=sd(mea)

# CBB bootstrap 
#   bl: block size
bsam_cbb=function(minta,bl)
{
  n=length(minta)
  bminta=c()
  minta2=c(minta,minta)
  bszam=ceiling(n/bl)
  hossz=1
  for (j in 1:bszam) {
    indul1=floor(runif(1,0,n))+1
    bminta = c(bminta,minta2[ (indul1:(indul1+bl-1)) ] )
    hossz=hossz+bl
  }
  bminta=bminta[1:n]
  return(bminta[bminta!=0])
}

# block bootstrap
n=length(yt)
mea2=rep(0,times=gen)         
k = b.star(yt,round=TRUE)[2]  # CBB, optimal block size
for (o in 1:gen) {
  bsample  = bsam_cbb(yt,k)
  mea2[o]=mean(bsample)           
}
becs2=sd(mea2)

# residual bootstrap
fitt=arima(yt,order=c(1,0,0))  # fits: p, d, q 
fitt
coef(fitt)
alfa_becs=coef(fitt)[1]
resids=fitt$residuals
mea3=rep(0,times=gen)         
for (o in 1:gen) {
  bsample = arima.sim(n=nnn, innov=sample(resids,replace=T), 
                       list(ar=c(alfa_becs)))
  mea3[o]=mean(bsample)           
}
becs3=sd(mea3)

# what is the true variance
valo=vector(length=nnn)
for (o in 1:nnn) {
  minta = arima.sim(n=nnn, list(ar=c(alfa)))
  valo[o]=mean(minta)           
}
valodi=sd(valo)
valodi2= sqrt(1/(nnn*(1-alfa^2))*(2*alfa^(nnn+1)-nnn*alfa^2-2*alfa+nnn )/
  (nnn*(alfa-1)^2))  # AR(1) formula, the exact "real" value

becs1-valodi2   # diff
becs2-valodi2
becs3-valodi2
100*(becs1-valodi2)/valodi2   # relative differences in %
100*(becs2-valodi2)/valodi2
100*(becs3-valodi2)/valodi2

# HW1: let us check the error for different b values. Which is the best?
#magyarul: ellenőrizzük a hibát különböző b értékekre. Melyik a legjobb?

# HW2: let us simulate from an AR(2) process (e.g. 0.5X_{t-1}+0.4X_t+\epsilon)
# and fit an AR(1) to it. What gives the parametric (residual-based) bootstrap based on this misspecified
#model?
#magyarul: szimuláljunk AR(2) folyamatból (pl. 0.5X_{t-1}+0.4X_t+\epsilon)
#és illesszünk rá AR(1) folyamatot. Mit ad a paraméteres (reziduálison alapuló bootstrap
#erre az esetre?

#MA
alfa=0.9
yt <- arima.sim(n=nnn, list(ma=c(alfa)))
par(mfrow=c(2,2))
plot(yt[1:100],type="l")
plot(ARMAacf(ma=c(alfa), lag=50), type="h")
plot(ARMAacf(ma=c(alfa), lag=50, pacf=TRUE), type="h")

b.star(yt,round=TRUE)  # optimal block size

# i.i.d. bootstrap
nn=gen
mea=rep(0,times=gen)
for (i in 1:nn)  mea[i]=mean(sample(yt,replace=T))
becs1=sd(mea)

# CBB bootstrap 
#   bl: block size
bsam_cbb=function(minta,bl)
{
  n=length(minta)
  bminta=c()
  minta2=c(minta,minta)
  bszam=ceiling(n/bl)
  hossz=1
  for (j in 1:bszam) {
    indul1=floor(runif(1,0,n))+1
    bminta = c(bminta,minta2[ (indul1:(indul1+bl-1)) ] )
    hossz=hossz+bl
  }
  bminta=bminta[1:n]
  return(bminta[bminta!=0])
}

# block bootstrap
n=length(yt)
mea2=rep(0,times=gen)         
k = b.star(yt,round=TRUE)[2]  # CBB, optimal block size
for (o in 1:gen) {
  bsample  = bsam_cbb(yt,k)
  mea2[o]=mean(bsample)           
}
becs2=sd(mea2)

# residual bootstrap
fitt=arima(yt,order=c(0,0,1))  # fits: p, d, q 
fitt
coef(fitt)
alfa_becs=coef(fitt)[1]
resids=fitt$residuals
mea3=rep(0,times=gen)         
for (o in 1:gen) {
  bsample = arima.sim(n=nnn, innov=sample(resids,replace=T), 
                       list(ma=c(alfa_becs)))
  mea3[o]=mean(bsample)           
}
becs3=sd(mea3)

# what is the true variance
valo=vector(length=nnn)
for (o in 1:nnn) {
  minta = arima.sim(n=nnn, list(ma=c(alfa)))
  valo[o]=mean(minta)           
}
valodi=sd(valo)
##########
#Task 4
##########

data = read.table("http://zempleni.elte.hu//reszv_ar.txt")
#library(np)

n=dim(data)[1]
par(mfrow=c(1,1))
plot(data[,1],type="l",ylim=range(data))
lines(data[,2],type="l",col=2)

##  a.)

acf(data[,1],lag.max=100)
acf(data[,2],lag.max=100)
n=dim(data)[1]
b1=b.star(data[,1],round=TRUE)[2]
b2=b.star(data[,2],round=TRUE)[2]
b1
b2

# not stationary, if there is a significant trend
magy=1:n

fit1=lm(data[,1]~magy)
summary(fit1)
plot(data[,1],type="l")
lines(predict(fit1),col=2)

fit2=lm(data[,2]~magy)
summary(fit2)
plot(data[,2],type="l")
lines(predict(fit2),col=2)  
# the linear trend is sign in both cases

# logreturns
logdata=cbind(diff(log(data[,1])),diff(log(data[,2])))
n=dim(logdata)[1]
magy2=1:n
par(mfrow=c(1,1))
plot(logdata[,1],type="l")
acf(logdata[,1],lag.max=100)
acf(logdata[,1]^2,lag.max=100)  # what does it show?
plot(logdata[,2],type="l")
acf(logdata[,2],lag.max=100)
acf(logdata[,2]^2,lag.max=100)  # what does it show?
fit1=lm(logdata[,1]~magy2)
summary(fit1)
fit2=lm(logdata[,2]~magy2)
summary(fit2)   # no sign trend
par(mfrow=c(2,1))
acf(logdata[,1],lag.max=100)
acf(logdata[,2],lag.max=100)
b1=b.star(logdata[,1],round=TRUE)[2]
b2=b.star(logdata[,2],round=TRUE)[2]
b1
b2

##  b.)

# block bootstrap to estimate the variance of the mean
gen=5000
mea.bl=array(dim=c(gen,2))
for (i in 1:gen) {
  mea.bl[i,1] = mean(bsam_cbb(logdata[,1],b1))
  mea.bl[i,2] = mean(bsam_cbb(logdata[,2],b2))         
}
sthiba=apply(mea.bl,2,sd)

##  c.)

bootatl=apply(mea.bl,2,mean)
atl=apply(logdata,2,mean)
hist(10^5*mea.bl[,1],breaks=30)
hist(10^5*mea.bl[,2],breaks=30) 
# diff. of the quantiles from the mean
alpha=0.05
10^5*(quantile(mea.bl[,1],c(alpha/2,1-alpha/2))-bootatl[1])
10^5*(quantile(mea.bl[,2],c(alpha/2,1-alpha/2))-bootatl[2])
   # bootstrap percentile interval

kvant=c(alpha/2,1-alpha/2)
sd1=quantile(mea.bl[,1]-atl[1],kvant)
confint1_v1=c(atl[1]-sd1[2],atl[1]-sd1[1])
sd2=quantile(mea.bl[,2]-atl[2],kvant)
confint2_v1=c(atl[2]-sd2[2],atl[2]-sd2[1])
10^5*confint1_v1
10^5*confint2_v1
# classic formula as bootstrap 
   # bootstrap-t interval, symmetric
seged=qnorm(kvant) # high sample size so a normal quantile can be used
confint1_v2=atl[1]+seged*sthiba[1]
confint2_v2=atl[2]+seged*sthiba[2]

10^5*confint1_v1   # just to see the values better
10^5*confint1_v2  
10^5*confint2_v1
10^5*confint2_v2