set.seed(12345)
sig2  = 1
tau2  = 0.25
ns    = c(seq(1000,10000,by=1000),seq(20000,100000,by=10000))
nn    = length(ns)
times = matrix(0,nn,2)
for (k in 1:nn){
  n    = ns[k]
  theta = rep(0,n)
  theta[1] = rnorm(1,0,sqrt(tau2))
  for (t in 2:n)
    theta[t] = rnorm(1,theta[t-1],sqrt(tau2))
  y = rnorm(n,theta,sqrt(sig2))

  #plot(y)
  #lines(theta,col=2)

  # Multivariate normal
#  time1 = system.time({
  Omega = matrix(1,n,n)
  for (i in 2:(n-1)){
    for (j in 1:i)
       Omega[i,j] = j
    for (j in (i+1):n)
       Omega[i,j] = i
  }
  Omega[n,]=1:n
  iOmega = solve(Omega)
  A = iOmega/tau2+diag(sig2,n)
  var = solve(iOmega/tau2+diag(sig2,n))
  mean = var%*%y
#})


# FFBS
time2 = system.time({
m0=0
C0=0
m = rep(0,n)
C = rep(0,n)
A    = (C0+tau2)/(C0+tau2+sig2)
m[1] = (1-A)*m0+A*y[1]
C[1] = A*sig2
for (t in 2:n){
  A    = (C[t-1]+tau2)/(C[t-1]+tau2+sig2)
  m[t] = (1-A)*m[t-1]+A*y[t] 
  C[t] = A*sig2
}
mn = rep(m[n],n)
Cn = rep(C[n],n)
for (t in (n-1):1){
  B = C[t]/(C[t]+tau2)
  Cn[t] = (1-B)*C[t]+B^2*Cn[t+1]
  mn[t] = (1-B)*m[t]+B*mn[t+1]
}
})


# par(mfrow=c(1,2))
# plot(y,ylim=range(y,theta,mean),xlab="Time",ylab="",pch=16)
# lines(mean,col=4,lwd=2)
# lines(mean-2*sqrt(diag(var)),col=4,lty=2,lwd=2)
# lines(mean+2*sqrt(diag(var)),col=4,lty=2,lwd=2)

# plot(y,ylim=range(y,theta,mean),xlab="Time",ylab="",pch=16)
# lines(mn,col=4,lwd=2)
# lines(mn-2*sqrt(Cn),col=4,lwd=2,lty=2)
# lines(mn+2*sqrt(Cn),col=4,lwd=2,lty=2)

# lines(m,col=2,lwd=2)
# lines(m-2*sqrt(C),col=2,lwd=2,lty=2)
# lines(m+2*sqrt(C),col=2,lwd=2,lty=2)

times[size,]=c(time1[1],time2[1])
}

par(mfrow=c(1,2))
plot(log10(sizes),times[,1],xlab="log10(n)",ylab="Time in seconds")
title("Matrix approach")
plot(log10(sizes),times[,1]/times[,2]/1000,xlab="log10(n)",ylab="Time rate matrix/ffbs (x1000)")
points(log10(sizes),times[,2],col=2)