------------------------------------------------
#  OLS-Schaetzung von Trend und Saison-Anteil  #
#            von Zeitreihendaten               #
------------------------------------------------

# viertel-jaehrliche BIP-Zahlen:
bipdata = read.table("C:/Users/Admin/desktop/BIPquarterly.csv",header=FALSE,sep=";")
bipdata
bip = bipdata[,2]
bip
summary(bip)
plot(bip)
plot(bip,type="l")      # hat offensichtlich eine 
                        # saisonale Komponente

# Bestimme die Saison- und Trend-Komponente 
# mit linearer Regression:
n = length(bip)    
n
t = 1:n                 # wie in BV4.1 indizieren wir die
                        # Zeitreihendaten mit ganzen Zahlen

# Regressoren fuer die Saison-Komponente: Da wir 
# quarterly data haben, ist die Periode T = 4:
X1 = cos(pi/2*t)
X2 = cos(pi*t)
X3 = sin(pi/2*t)

X = cbind(X1,X2,X3)

# Regressoren fuer die Trend-Komponente: 
# einfache Polynome bis zur Ordnung m:
m = 8
for(k in 1:m)
{
  X = cbind(X,t^k)
}

# damit ist alles upgesetted und wir koennen 
# die lineare Regression durchfuehren:

res = lm(bip ~ X)

plot( bip , type="l" )
lines( res$fit , col=2 )

# let's isolate the trend and season-component:
beta = res$coef
beta

trendcomp = beta[1] * rep(1,n)
for(k in 1:m)
{
  trendcomp = trendcomp + beta[k+4]*t^k
}

# let's have a look:
lines( trendcomp , col=3 )

seasoncomp = beta[2]*X1 + beta[3]*X2 + beta[4]*X3
plot( seasoncomp , type="l" )

# let's check the regression fit:
TestRegFit = trendcomp + seasoncomp
plot( bip , type="l" , ylim=c(-50,650) )
lines( res$fit , col=2 )
lines( TestRegFit , col=3 )
lines( trendcomp , col=2 )
lines( seasoncomp ,col=4 )


# the actual procedure in BV4.1 is more subtle, 
# but roughly the pictures are the same: 
# let's have a look  ->  BV4.1