################################################# ### Numerische Optimierung ################################################# # In VWL, ?konometrie und Statistik ist h?ufig Maximierung # oder Minimierung von Funktionen gefordert. Bei einigen Funktionen # (statistische Zielfunktionen wie Likelihood oder Residuenquadratsumme # nichtlinearer Regressionen) ist das Optimierungsproblem analytisch # nicht l?sbar, es existiert keine (algebraische) Formel f¨¹r den Sch?tzer (wie etwa bei OLS) # Ausweg: Es ist m?glich, die Funktion NUMERISCH bzw. allgemeiner ITERATIV (in mehreren Schritten) # zu optimieren. Dazu muss i.A. nur der Funktionswert evaluierbar sein. # Achtung: Bei numerischen Algorithmen muss auf Eindeutigkeit der Nullstelle, des Optimums, # des Fixpunkts usw geachtet werden. U. a. k?nnen unterschiedliche Startwerte zu # unterschiedlichen Ergebnissen f¨¹hren (s. unten) # Funktionen optimize (univariat), nlm() und optim() (multivariat) # -> siehe Fortgeschrittene ?konometrtie f¨¹r die Darstellung der Methoden # Hier nur optim verwendet (am allgemeinsten!) ?optim #### Beispiel 1: Einfache univariate Funktion f <- function(x) 0.2*(x-1)^4-0.7*x^2-0.3*x+5 x <- seq(0,4,0.01) plot(x,f(x),type="l") f(2) # Minimum bestimmen # par = Startwert, fn=functionname der Zielfunktion, lower = Untergrenze bei Brent, # upper = Obergrenze ?optim optim(par=NA, fn=f, lower=-10, upper=10, method="Brent") abline(v = 2.726503) ### Aufgabe H1 # Modifizieren Sie die Funktion f, so dass vor der Ausgabe des Funktions- # werts mit der print-Funktion der x-Wert und der Funktionswert ausgegeben # wird. F¨¹hren Sie nun die Optimierung durch. Wie erkl?ren Sie sich das? ### Aufgabe H2 f(x) = exp(x) + log(x) # Lesen Sie sich die Wikipediaeinleitung zur Bisektion durch. # Schreiben Sie nun einen Algorithmus, der Nullstellen einer Funktion auf einem Inter- # vall [-4,6] sucht und testen Sie ihn an der Funktion 3*x-4 # # Hinweis: # Schreiben Sie eine while-Schleife mit der Abbruchbedingung, dass der # Unterschied zweier Iterationsschritte kleiner 0.01 ist. # Berechnen Sie in der Schleife den Mittelpunkt des Intervalls, # vergleichen Sie den Wert des Mittelpunkts mit den # beiden anderen Werten der Intervallgrenzen, w?hlen Sie nun f¨¹r den n?chsten Schritt das Intervall, # bei dem sich die Nullstelle befinden muss -> Vorzeichenwechsel! #### Wieso Optimierung in der ?konometrie? # Viele Sch?tzer basieren auf der Minimierung einer Funktion von den Residuen # Beispiel: Nichtlineare Regression (siehe Fortgeschrittene ?k.) ## Regressionsmodell mit zwei Parametern a1 und a2 a <- c( 1 , 0.6 ) n <- 200 x <- rchisq(n,df=3) u <- rnorm(n,sd=0.5) # Nichtlineares Modell y = a x^b+u y <- a[1]*x^a[2] + u # Meist verwendet man stattdessen log(y) = c0 + c1 log(x) + e DATA <- data.frame(y,x) plot(DATA$x,DATA$y) ## Residuenquadratsumme zu minimieren RSS <- function(a,DATA){ u_tilde <- DATA$y - a[1]*DATA$x^a[2] return(sum(u_tilde^2)) } ## Zielfunktion anschauen f¨¹r verschiedene Werte von a length.grid <- 50 a1_grid <- seq(-1,3,length.out=length.grid) a2_grid <- seq(0,1.5,length.out=length.grid) a_grid <- as.matrix(expand.grid(a1=a1_grid,a2=a2_grid)) RSS_grid <- apply(a_grid,1,RSS,DATA=DATA) RSS_grid <- matrix(RSS_grid,length.grid,length.grid) # Nicht viel zu sehen contour(a1_grid,a2_grid,RSS_grid,xlab="a1",ylab="a2",nlevels=100) # -> Logarithmieren contour(a1_grid,a2_grid,log(RSS_grid),xlab="a1",ylab="a2",nlevels=100) # 3D-Plot ?persp persp(RSS_grid,xlab="a1",ylab="a2",phi=10,theta=-40,ticktype="detailed") # damit das Minimum sp?ter auch sinn macht! persp(log(RSS_grid),xlab="a1",ylab="a2",phi=10,theta=-40,ticktype="detailed") ## Numerische Minimierung # par = Startwert # fn = Zielfunktion est <- optim(par=c(1,1), fn=RSS, DATA=DATA) points(x = est$par[1], y = est$par[2] , col = "red" , pch = 20) # Fazit: Am Tripel (0.99, 0.61, 45.14) = (a1, a2, f(a1,a2)) # liegt das Minimum der Residuenquadratsumme, da am Punkt (a1, a2) # die Funktion RSS das Minimum, den gesch?tzten Wert 45.14 annimmt # Allgemein: # Hat man eine Funktion f: R^m -> R^n, so hat der Graph Dimension n+m # In diesem Fall: R^m=R^2 und R^n=1, also dreidimensional ?optim ## Bei nichtlinearer Regression einfacher: nls() est_nls <- nls(y~a1*x^a2, start=c(a1=1,a2=1),data=DATA) summary(est_nls) #### Beispiel 3: Startwerte sind wichtig! # Ziel: Maximieren der Funktion f f <- function(x) (2*cos(x[1]/0.5)-0.3*x[1]^2)+cos(x[2]^2)-0.5*(x[2]-2)^2 # max statt min: control-Argument benutzen (oder g=-f minimieren) optim(par=c(1,1),f,control=list(fnscale=-1)) optim(par=c(-0.5,3),f,control=list(fnscale=-1)) # Andere Startwerte... # Funktion betrachten l.grid <- 50 grid.x <- as.matrix(expand.grid(x1=seq(-5,5,length.out=l.grid),x2=seq(-5,5,length.out=l.grid))) grid.f <-apply(grid.x,1,f) matrix.grid.f <- matrix(grid.f ,l.grid) image(matrix.grid.f,xlab="x1",ylab="x2") contour(matrix.grid.f,nlevels=100,add=TRUE) persp(x = seq(-5,5,length.out=l.grid), y = seq(-5,5,length.out=l.grid), z = matrix.grid.f, xlab = "x1", ylab = "x2", theta = 60, ticktype = "detailed") # Und Startwerte geschickt w?hlen (x.start <- grid.x[which.max(grid.f),]) optim(par=x.start,f,control=list(fnscale=-1)) ### Aufgabe H3 # Maximieren Sie die Funktion f f <- function(x) (2*cos(x[1]/0.5)-0.3*x[1]^2)+cos(x[2]^2)-0.5*(x[2]-2)^2 # Gehen Sie davon aus, dass das Maximum (x1, x2) -5 Var(x[t]|I[t-1]) = h[t]^2 Bedingte Varianz ?ndert sich (ist von Vorperioden abh?ngig:) # h[t]^2 = c + g[1]*h[t-1]^2 + ... + g[q]*h[t-q]^2 + a[1]*x[t-1]^2 + ... + a[p]*x[t-p]^2 ### Definiere die Log Likelihood Funktion ### theta = Parametervektor ### x = Vektor (univariate Zeitreihe) LogLik.GARCH <- function(theta, x) { c <- theta[1] # Konstante in GARCH Gleichung a <- theta[2] # ARCH Parameter g <- theta[3] # GARCH Parameter # Unzul?ssige Parameterwerte ausschlie?en (-> Inf als Wert) if (sum(a)+sum(g)>1) return(NA) # Instabil -> nicht zugelassen if (any(theta<0)) return(NA) # Nichtnegativit?t verletzt -> nicht zugelassen n <- length(x) # Stichprobengr??e h_sq <- rep(var(x),n) # Vektor der bedingten Varianzen (am Anfang: unbedingte Var.) ll <- numeric(n) # Vektor der (MINUS) Log Likelihoods pro Beobachtung # Iterativ Bedingte Varianz und minus Likelihood berechnen for (i in 2:n) { h_sq[i] <- c+a*x[i-1]^2+ g*h_sq[i-1] ll[i] <- -0.5*(log(2*pi)+log(h_sq[i])+x[i]^2/h_sq[i]) # minus Likelihood } # Funktionswert = Log Likelihood return(sum(ll)) } ### Aufgabe H4 # a) Suchen und finden Sie eine m?glichst lange Zeitreihe zu t?glichen Preisen # von Getreidefutures. # b) Finden Sie f¨¹r x='Vektor der log-Preis?nderungen' (Renditen) # den Parametervektor, der die oben definierte Funktion maximiert # -> Maximum Likelihood Sch?tzer # c) Berechnen Sie die bedingten Varianzen (h_sq) # d) Plotten Sie untereinander Renditen und bedingte Varianzen # e) Prognostizieren (berechnen) Sie die HEUTIGE bedingte Varianz der Renditen! # Heute sei Aug 2012, letztes vorhandener Wert von Mar 2012 -> 5-step forecast # f) Ist diese gr??er/kleiner als "¨¹blich"? ### Aufgabe H5 # Im File FunktMIN.RData sind die Vektoren y, x und z gespeichert. # Betrachten Sie das Regressionsmodell # y = beta[1]+beta[2]*x+beta[3]*z^beta[4] + u # a) Geben Sie eine Funktion an, die aus den Daten und beliebigen Werten von beta # die Residuenquadratsumme berechnet. # b) Berechnen Sie den Wert von beta, der diese SSR minimiert. # Startwerte aus linearen Regression (mit beta[4]=1) # c) Versuchen Sie, die Funktion f im File FunktMIN.RData # zu minimieren. Geben Sie auch die Hessematrix an.