# The Additive Model for a Time Series (Solution to Homework)

 
# Ex. 1.1 The Mitscherlich function  with different values of the parameters
t <- seq(0,6,length=100)
Mfun <- function(t,b1,b2,b3) b1 + b2*exp(t*b3)
plot(range(t),c(0,2),type="n",xlab="t",ylab="Mfun")
lines(t,Mfun(t,1,-1,-0.5))
lines(t,Mfun(t,1,1,-0.5),col="red")
lines(t,Mfun(t,1,-1,-1),col="blue")
lines(t,Mfun(t,1,1,-1),col="green")

# Ex. 1.5 Fitting the logistic model to 
population2 <- read.table("population2.txt",col.names=c("year","t","pop"))
population2
dim(population2)
# Regression fit
z.t <- 1/population2$pop[-10]
z.p <- 1/population2$pop[-1]
plot(z.t~z.p)
x <- 1:10
y <- population2$pop
fit <- lm(z.t~z.p)
a <- coef(fit)[1]; b <- coef(fit)[2]
b1 <- log(b)
b3 <- (1-exp(-b1))/a
b2 <- sum(((b3-y)/y)*exp(-b1*x))/sum(-2*b1*x)
b1;b2;b3
# nls fit
fit.pop1 <- nls(pop ~ b3/(1+b2*exp(-t*b1)),data = population2,
                 start = list(b1=b1,b2=b2,b3=17))
fit.pop1
# Comparison
# The linear fit does not make sense. b3 is negative since the intercept a is.


# Ex 1.6 (Income Data). The (accumulated) annual average 
# increases of gross and net incomes in thousands DM (deutsche mark) 
# in Germany, starting in 1960.
income <- read.table("income.txt")
colnames(income) <- c("year","t","gross","net")
income
plot(net~t,data=income)
y <- income$net[-1]
x <- income$t[-1]
# Regression fit
fit.reg <- lm(log(y) ~ log(x))
summary(fit.reg)
coef(fit.reg)
b1 <- coef(fit.reg)[2]
b2 <- exp(coef(fit.reg)[1])
b1
b2
y.hat <- b2*x^b1
1-sum((y-y.hat)^2)/sum((y-mean(y))^2)
# nls fit
fit.nls <- nls(y ~ b2*x^b1,start = list(b1=b1,b2=b2))
1-sum((y-predict(fit.nls))^2)/sum((y-mean(y))^2)
# Comparison
# By definition the nls fit is better, but the regression  fit is not bad.


# Ex. 1.7
Still to come