# Init by deleting all variables and functions
rm(list=ls())


#----------------------------------------------------------------
# Exercise on RLS

# Generate some data from a linear model
n <- 200
x <- runif(n)
beta0 <- 2
beta1 <- -3
y <- beta0 + beta1 * x + rnorm(n, sd=0.1)

# Try to estimate the parameters
fit <- lm(y ~ x)
plot(x, y)
abline(fit)
fit

# Change the coefficients and generate more data
x1 <- runif(n)
beta0 <- -2
beta1 <- 3
y1 <- beta0 + beta1 * x1 + rnorm(n, sd=0.1)

## Put together, such that at the start of the the second part, i=201 the coefficients
## change and have different values
y <- c(y, y1)
x <- c(x, x1)

# See clearly there are two "regimes"
plot(x, y)

# Here you only see the change in mean
plot(y)

# Fit a linear regression on the binded data
lm(y ~ x)

## A function for fitting a recursive least squares estimation
rls <- function(lambda){
  ## Single inupt RLS
  ## Must have y, X, lambda
  ## Number of obs
  n <- nrow(X)
  ## Number of parameters
  p <- ncol(X)
  ## For keeping the parameters
  Thetahat <- matrix(NA, ncol=p, nrow=n)
  ## For keeping predictions
  yhat <- rep(NA, n)
  ##---- The recursive variables
  theta <- matrix(rep(0,p), ncol = 1)
  ## Start value for the parameter covariance P
  P <- diag(10000, p)
  ## Or the inverse in the RLS
  R <- solve(P)
  ## Iterate through and estimate the parameters
  for(i in 1:n){
    xi <- matrix(X[i, ])
    yhat[i] <- t(xi) %*% theta
    ## Update the parameter estimate
    R <- lambda * R + xi %*% t(xi)
    theta <- theta + solve(R, xi) %*% (y[i] - t(xi) %*% theta)
    ## Keep output
    Thetahat[i, ] <- t(theta)
  }
  ## Return a list
  L <- list()
  L$yhat <- yhat
  L$Thetahat <- Thetahat
  return(L)
}

## Fit a Recursive Least Squares (RLS)
## Forgetting factor
lambda <- 0.99
## Design matrix
X <- cbind(matrix(1, ncol=1, nrow=length(x)), matrix(x, ncol=1))
## Run it
val <- rls(lambda)
# Plot the estimated parameters
colnames(val$Thetahat) <- c("beta0","beta1")
plot.ts(val$Theta)

obj <- function(lambda){
  val <- rls(lambda)
  sum((val$yhat - y)^2)
}
obj(0.1)
obj(0.99)
optimize(obj, c(0.01,1))

##----------------------------------------------------------------