install.packages(c("dynlm","vars"),repos="http://cran.r-project.org")

#################################################
### Basic
#################################################

# Aus empirischer Verteilungsfunktion der Daten (Residuen) wird die Verteilung
# von Teststatistiken oder anderer abgeleiteter Gr??en simuliert.
# Idee: Annahme, die Stichprobe entspreche der ganzen Population und nutze nun 
# diese Population um durch Ziehen mit Zur?cklegen neue Stichproben zu erzeugen. 
# Siehe Mehoden der ?konometrie, Davidson/McKinnon (2004, Abschn 4.6)


# Simuliere AR(1) Prozess (mit nichtnormal verteiltem Fehler)
# y(t) = a*y(t-1) + u(t)

set.seed(400)
a <- 0.9
n <- 150
u <- rt(n,df=5)
?filter
y <- filter(u,a,"recursive")
plot.ts(y)

# Sch?tze a
install.packages("dynlm")
library(dynlm)
est.ar <- dynlm(y~L(y,1))


(a_hat <- est.ar$coef[2])
(u_hat <- est.ar$resid)

### Test auf H_0: a = 0.9

# t-Test (berechne p-value), hier nur asymptotisch g?ltig
t_stat <- (summary(est.ar)$coef[2,1]-0.9)/summary(est.ar)$coef[2,2]
pnorm(abs(t_stat),lower.tail=FALSE)*2


# (Manchmal bessere) Alternative: Bootstrap-Test
# Ziehe aus Residuen (mit zur?cklegen)
# Generiere Reihen mit H_0 Wert
# Berechne jeweils t-Statistik

Reps <- 999
a_sim <- numeric(Reps)
t_sim <- numeric(Reps)

for (i in 1:Reps){
	u_sim <- sample(u_hat,n,replace=TRUE)     # Ziehen aus den Residuen (mit Zur?cklegen!)
	y_sim <- filter(u_sim,0.9,"recursive")    # Nullhypothesenwert verwenden!
	est_sim <- dynlm(y_sim~L(y_sim,1))
	a_sim[i] <- est_sim$coef[2]
	t_sim[i] <- (a_sim[i]-0.9)/summary(est_sim)$coef[2,2]
}

# Berechne P-Value
mean(abs(t_sim)>abs(t_stat))

plot(density(t_sim),main="Asymptotic and bootstrap distribution of t",xlab="t-Statistic")
lines(seq(-4,4,0.1),dnorm(seq(-4,4,0.1)),col="orange")
legend("topright",c("Bootstrap","Asymptotic"),lty=c(1,1),col=c(1,"orange"))


# Konfidenzintervalle f?r Koeffizienten  (L?tkepohl, 2005; Appendix)
# der MA-Darstellung (=Impuls-Antwort-Funktion)
# y(t) = u(t) + phi(1) u(1) +  phi(2) u(2) + ...

# Wahre Impulsantworten (IRs)
(ir_true <- a^(0:30))
plot(0:30,ir_true,type="l",ylab="")

# Gesch?tzte IRs
ir_est <- a_hat^(0:30)
lines(0:30,ir_est,lty=2,lwd=2)

# Berechne simulierte IRs
# Plotte gesch?tzte und einige simulierte IRs
# Resampling
for (i in 1:Reps)
{
	u_sim <- sample(u_hat,n,replace=TRUE)
	y_sim <- filter(u_sim,a_hat,"recursive")     # Hier mit gesch?tztem Wert
	est_sim <- dynlm(y_sim~L(y_sim,1))
	a_sim[i] <- est_sim$coef[2]
}

IR_sim <- matrix(NA,Reps,31)
for (j in 1:Reps) IR_sim[j,] <- a_sim[j]^(0:30)

plot(c(0,30),c(0,1),type="n")
for (i in sample(1:Reps,90)) lines(0:30,IR_sim[i,],col=i,lty=3)
lines(0:30,ir_est,type="l",lty=1,lwd=2)
legend("topright","Gesch?tzte IAF",lty=1,lwd=2)


# Berechne Konfidenzintervalle (Perzentile der simulierten Verteilung)

ir_lower <- apply(t(IR_sim),1,quantile,probs=0.025)
ir_upper <- apply(t(IR_sim),1,quantile,probs=0.975)

plot(0:30,ir_est,ylim=c(0,1.02),type="l",
	xlab="Horizon",ylab="Impulse Response",
	main="Estimated Impulse Responses with 95% Confidence Intervals")
lines(0:30,ir_lower,lty=2)
lines(0:30,ir_upper,lty=2)



#################################################
### Anwendung Bootstrap-Konfidenzintervalle:
### "Technology shocks and hours worked"
#################################################

setwd("g:/R Kurs/Daten2")
HOABS <- read.table("HOABS.txt",skip=11,header=TRUE)
OPHPBS <- read.table("OPHPBS.txt",skip=17,header=TRUE)
CNP16OV <- read.table("CNP16OV.txt",skip=17,header=TRUE)

DATA <- merge(HOABS, OPHPBS, by="DATE")
DATA <- merge(DATA, CNP16OV, by="DATE")
names(DATA) <- c("DATE","HOABS","OPHPBS","CNP16OV")
DATA$productivity <- log(DATA$OPHPBS)
DATA$d.productivity <- c(NA,diff(log(DATA$OPHPBS)))
DATA$hours <- log(DATA$HOABS/DATA$CNP16OV)

head(DATA)
plot.ts(DATA[,5:7], main="Productivity and Hours Worked")
y <- as.matrix(na.omit(DATA[,6:7]))

# Frage: F?hrt ein positiver Technologieschock in der Folge zu mehr
#        geleisteten Arbeitsstunden (RBC) oder zu weniger
#        (technologisch bedingte Arbeitslosigkeit)?

# -> Impuls-Antwort Analyse in VARs

# Vektor Autoregressives Modell (vgl. Quantitative WiFo; L?tkepohl, 2005)
# -> Regression einer multivariaten Zeitreihe auf eigene verz?gerte Werte
# (soviele Regressionsgleichungen wie Anzahl der Zeitreihen)
# y[t] = A1 y[t-1] + A2 y[t-2] + ... + Ap y[t-p] + u[t]

# Sch?tzung der VAR-Koeffizienten A1, A2, ... mit OLS:
p <- 4
Y.lags <- embed(y,p+1)
Y.lags[1:10,]
est_ols <- lm(Y.lags[,1:2]~Y.lags[,3:((p+1)*2)])
summary(est_ols)
B <- est_ols$coef[-1,]
A <- array(t(B),c(2,2,p))

# Berechnung der (Vektor-) Moving Average Parameter Mj (= Impuls-Antwort-Funktion)
# aus
# y[t] = u[t] + M1 u[t-1] + M2 u[t-2] + ...
h <- 100
AA <- rbind(matrix(A,2),cbind(diag(2*(p-1)),matrix(0,2*(p-1),2)))
MM <- array(0,c(p*2,p*2,h))
MM[,,1] <- diag(p*2)
for (i in 2:h) MM[,,i] <- MM[,,i-1]%*%AA
M <- array(MM[1:2,1:2,],c(2,2,h))

### Plot der Impuls-Antwort Funktion (Reaktion der Arbeitsstunden auf Schocks der
#   Arbeitsproduktivit?t)
plot(M[2,1,],type="l",xlab="Horizon",ylab="Reaction")

### Aufgabe:
### Geben Sie Bootstrap-Konfidenzintervalle f?r diese Impuls-Antwortfunktion an
# Hinweis: Ein VAR-Prozess l?sst sich wie folgt aus Innovationen u simulieren:
# for (j in (p+1):(n+p))
#      y[j,] <- A[,,1]%*%y[j-1,]+A[,,2]%*%y[j-2,]+A[,,3]%*%y[j-3,]+A[,,4]%*%y[j-4,]+u[j,]


### Probleme:

### Problem 1: Fehler aus der reduzierten Form, u[t], stellen keine "eigenst?ndigen"
### bzw. strukturellen Schocks auf Produktivit?t bzw. Arbeitsstunden dar ->
### Strukturelle VAR Modelle

### Problem 2: Verbesserte Methoden zur Berechnung von Bootstrap-Konfidenzintervallen
### vgl. L?tkepohl (2005); etwa Hall-Methode

### Problem 3: Sind Arbeitsstunden station?rer Prozess? Oder Unit Root?
### -> Laufende Debatte in der Literatur...

### Strukturelles VAR a la Blanchard&Quah: (Strukturelle) Arbeitszeitschocks haben keine langfristige
### Wirkung auf Produktivit?t
library(vars)
redform <- VAR(y,ic="AIC",lag.max=8)
strucform <- BQ(redform)
ir_struc <- irf(strucform,n.ahead=100,runs=1000)
str(ir_struc)

plot(ir_struc$irf$d.productivity[,"hours"],
  ylim=c(-0.001,0.018),type="l",
  main="Reaction of hours worked to STRUCTURAL productivity shock",
  ylab="Responses",xlab="horizon")
lines(ir_struc$Upper$d.productivity[,"hours"],lty=2)
lines(ir_struc$Lower$d.productivity[,"hours"],lty=2)
abline(h=0,lty=3)



#################################################
### Meta-Simulation
#################################################

# Durchf?hren einer Simulation, bei der Eigenschaften einer
# Simulationsmethode ?berpr?ft werden.
# Etwa: Hat der Bootstrap bessere Eigenschaften als asymptotische Tests?


### Aufgaben I1
# Wiederholen Sie die Bootstrap-Prozedur im AR(1) (oder einem anderen) Beispiel
# F?hren Sie jeweils einen Bootstrap- und einen asymptotischen (t-) Test bzgl.
# a durch. Wie gut ist der Bootstrap bei wahrer Nullhypothese?



### Aufgaben I2
# Plotten Sie die Ablehnungsh?ufigkeiten des Bootstrap- und asymptotischen
# Tests in Abh?ngigkeit von a.


### N?tzliche Funktion, die aus den Daten die Impuls-Antwort-Funktion sch?tzt

IR_estimate <- function(y,p,h){
  Y.lags <- embed(y,p+1)
  A <- array(t(lm(Y.lags[,1:2]~Y.lags[,3:((p+1)*2)])$coef[-1,]),c(2,2,p))
  AA <- rbind(matrix(A,2),cbind(diag(2*(p-1)),matrix(0,2*(p-1),2)))
  MM <- array(0,c(p*2,p*2,h))
  MM[,,1] <- diag(p*2)
  for (i in 2:h) MM[,,i] <- MM[,,i-1]%*%AA
  M <- array(MM[1:2,1:2,],c(2,2,h))
  return(M[2,1,])
}

### Berechne Residuen aus OLS-Sch?tzung
u <- est_ols$resid
n <- nrow(u)

### Sch?tze Impuls-Antwort-Funktion aus simulierten Daten
B <- 1000
IR_sim <- matrix(0,h,B)

for (i in 1:B){
  us <- u[sample(1:n,n+p,replace=TRUE),] # Ziehe aus Residuen mit Zur?cklegen
  ys <- us
  for (j in (p+1):(n+p))                # Rekursiv VAR Prozess simulieren
    ys[j,] <- A[,,1]%*%ys[j-1,]+A[,,2]%*%ys[j-2,]+A[,,3]%*%ys[j-3,]+A[,,4]%*%ys[j-4,]+us[j,]
  # Jeweils Impuls-Antwort-Funktion sch?tzen und speichern
  IR_sim[,i] <- IR_estimate(y=ys,p=4,h=h)
}

### Plot der simulierten IAFen
for (i in 1:B) lines(IR_sim[,i],lwd=0.4,lty=2)
lines(M[2,1,],col=2,lwd=2)

### Konfidenzintervalle berechnen
ir_lower <- apply(IR_sim,1,quantile,probs=0.025)
ir_upper <- apply(IR_sim,1,quantile,probs=0.975)

plot(0:(h-1),M[2,1,],ylim=c(0,1.2),type="l",
     xlab="Horizon",ylab="Reaction of working hours to technology shock",
     main="Estimated impulse responses with 95% confidence intervals")
lines(0:(h-1),ir_lower,lty=2)
lines(0:(h-1),ir_upper,lty=2)