Nothing
#' VBV.decomposition - decompose a time series with VBV
#'
#' VBV.decomposition - decompose a time series with VBV
#'
#' @param n number of observation points. Internally this will be transformed to
#' seq((-(n-1)/2, (n-1)/2, 1)
#'@param p maximum exponent in polynomial for trend
#'@param q.vec vector containing frequencies to use for seasonal component, given as integers, i.e. c(1, 3, 5) for 1/2*pi, 3/2*pi, 5/2*pi (times length of base period)
#'@param grundperiode base period in number of observations, i.e. 12 for monthly data with yearly oscillations
#'@param lambda1 penalty weight for smoothness of trend
#'@param lambda2 penalty weight for smoothness of seasonal component
#' lambda1 == lambda2 == Inf result in estimations of the original Berliner Verfahren
#' @return list with the following components:
#' \item{trend}{A function which returns the appropriate weights if applied to a point in time}
#' \item{saison}{A function which returns the appropriate weights if applied to a point in time}
#' \item{A, G1, G2}{Some matrices that allow to calclate SSE etc. Exposed only reuse their calculation. See the referenced paper for details.}
VBV.decomposition <-
function(n, p , q.vec, grundperiode, lambda1, lambda2)
{
### Die Funktion VBV.zerlegung implementiert das verallgemeinerte
### Berliner Verfahren auf Grundlage des Berliner Verfahrens in
### der Matrixdarstellung.
### lambda1 = lambda2 = Inf ergeben die Schäter des Berliner Verfahrens
### in der Grundversion
### n: Länge des Zeitvektors, wird umgewandelt in seq((-(n-1)/2, (n-1)/2, 1)
### p: maximaler Grad des Trendpolynoms
### q.vec: Vektor der zu betrachtenden Schwingungen als ganze Zahlen
### z.B. c(1, 3, 5) für 1/2*pi, 3/2*pi, 5/2*pi (*Grundperiodenlänge)
### grundperiode: Grundperiode in Anzahl Beobachtungen. Z.B. 12 bei
### Monatsdaten mit jährlichen Schwingungen
### lambda1: Strafgewicht für den Trend
### lambda2: Strafgewicht für die Saison
### lambda1 = lambda2 = Inf ergeben das Berliner Verfahren!
### Rückgabewert ist eine Liste mit den Komponenten:
### trend: Funktion(!) die angewendet auf einen Zeitpunkt aus dem Zeitbereich den
### Gewichtsvektor liefert
### saison: Funktion(!) die angewendet auf einen Zeitpunkt aus dem Zeitbereich den
### Gewichtsvektor liefert
### A, G1, G2: Matrizen, die erlauben Fehlerquadrate etc zu berechnen.
t.vec <- seq(-(n-1)/2, (n-1)/2, 1)
f10 <- function(t, p, q){
c(t^(0:(p-1)), rep(0,2*q) )
}
f02 <- function(t, p, q.vec, grundperiode){
## Seite 43
q <- length(q.vec)
erg <- rep(0,p+2*q)
erg[seq(p+1,p+2*q-1,2)] <- cos(2*pi*t*q.vec/grundperiode)
erg[seq(p+2,p+2*q,2)] <- sin(2*pi*t*q.vec/grundperiode)
erg
}
g1.plus <- function(tk,t,p)
{
if (t <= tk) return (0)
(t-tk)^(2*p -1)
}
g1 <- function(t, t.vec, p)
{
sapply(t.vec, g1.plus, t=t, p=p)
}
g2.plus <- function(tk,t,grundperiode,q.vec)
{
if (t <= tk) return(0)
r.vec <- q.vec/mean(q.vec)
r2.vec <- r.vec^2
d.vec <- 1/r.vec
c.vec <- 1/r.vec
q <- length(q.vec)
for (i in 1:q)
{
d.vec[i] <- d.vec[i] - 4*r.vec[i] * sum(1/(r2.vec[-i]-r2.vec[i]))
c.vec[i] <- c.vec[i] / prod(r2.vec[-i]-r2.vec[i])
}
sum(c.vec^2*( d.vec * sin( 2*pi* (t-tk)*q.vec/grundperiode) -
mean(q.vec) * 2 * pi * (t -tk)/grundperiode *
cos(2*pi*(t-tk)*q.vec/grundperiode) )
)
}
g2 <- function(t, t.vec, q.vec, grundperiode)
{
sapply(t.vec, g2.plus, t=t, q.vec=q.vec, grundperiode=grundperiode)
}
FGES <- t(
sapply(t.vec, f10, p=p, q=length(q.vec)) +
sapply(t.vec, f02, p=p, q.vec=q.vec, grundperiode=grundperiode)
)
Bstar <- qr.solve( (t(FGES)%*%FGES), tol=1e-17) %*%t(FGES)
Astar <- diag(1, length(t.vec)) - FGES%*%Bstar
G1 <- sapply(t.vec, g1, t.vec=t.vec, p=p)
G2 <- sapply(t.vec, g2, t.vec=t.vec, q.vec=q.vec, grundperiode=grundperiode)
GGES <- -t( G1/lambda1 + G2/lambda2 )
A <- qr.solve(diag(1, n) - Astar%*%GGES, tol=1e-17)%*%Astar
## A1 <- qr.solve(diag(1, n) - Astar%*%(diag(1,n)-G1), tol=1e-17)%*%Astar
## A2 <- qr.solve(diag(1, n) - Astar%*%(diag(1,n)-G2), tol=1e-17)%*%Astar
B <- Bstar%*% (diag(1, n)+ GGES %*% A)
list(
trend = function(t) {
f10(t,p,length(q.vec))%*% B +
g1(t, t.vec,p)%*%((1/lambda1)*A)
},
saison = function(t) {
f02(t,p,q.vec,grundperiode) %*% B +
g2(t, t.vec, q.vec, grundperiode) %*% ((1/lambda2)*A)
},
A = A , Astar = Astar,
G1 = G1, # sapply(t.vec, g1, t.vec=t.vec, p=p),
G2 = G2 #sapply(t.vec, g2, t.vec=t.vec, q.vec=q.vec, grundperiode=grundperiode)
)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.