tvec<-seq(0.1,2,by=0.1)
ey<-1+0.1*tvec+0.2*log(tvec)
plot(tvec,ey,type="l")

n<-200
tv<-seq(0.2,2.0,by=0.2)
Tv<-length(tv)

Xi<-cbind(rep(1,Tv),tv,log(tv))
be<-c(1,0.1,0.2)
p<-3

set.seed(33)
sigu<-0.04
sige<-0.01
u<-rnorm(n)*sigu
ID<-numeric()
Y<-numeric()
X<-numeric()
pdf("GrowthData.pdf",paper="USr",width=11,height=9)
plot(tvec,ey,type="n",xlab="t",ylab="y",ylim=range(0.5,1.5),xlim=range(0,2))
for(i in 1:n){

	Yi<-Xi %*% be + u[i] + rnorm(Tv)*sige
	lines(tv,Yi)

	ID<-c(ID,rep(i,Tv))
	X<-rbind(X,Xi)
	Y<-c(Y,Yi)

}
dev.off()

######################################################
#Establish well calibrated priors
sigma.0<-10
a.u<-20
b.u<-1
hist(1/rgamma(10000,a.u,b.u));box()

a.e<-5
b.e<-0.05
hist(1/rgamma(10000,a.e,b.e));box()

######################################################
library(MCMCpack)
library(mnormt)
nburn<-100
nthin<-20
nsamp<-2000
nits<-nburn+nthin*nsamp

usd<-0.02

old.u<-rep(0,n)
old.uvec<-rep(old.u,each=Tv)
old.beta<-as.numeric(coef(lm(Y~X-1)))
old.muvec<-X %*% old.beta + old.uvec
old.res<-Y-old.muvec
old.ssq<-(old.res)^2


old.sigma.e<-0.02
old.sigma.u<-0.05
old.clike<--0.5*n*Tv*log(old.sigma.e^2)+sum(-0.5*(Y-old.muvec)^2)/old.sigma.e^2
old.ilike<-apply(matrix(-0.5*old.ssq/old.sigma.e^2,nrow=n,byrow=T),1,sum)

old.uilike<-dnorm(old.u,0,old.sigma.u,log=T)
old.ulike<-sum(old.uilike[-1])

old.prior.beta<--0.5*p*log(old.sigma.e^2*sigma.0^2)-0.5*sum(old.beta)^2/(old.sigma.e^2*sigma.0^2)

old.prior.sigu<-log(dinvgamma(old.sigma.u^2,a.u,b.u))

ico<-0
growth.beta<-matrix(0,nrow=nsamp,ncol=p)
growth.sigmas<-matrix(0,nrow=nsamp,ncol=2)
growth.us<-matrix(0,nrow=nsamp,ncol=n)

t0<-proc.time()[3]
for(iter in 1:nits){

	#Update the us using Metropolis
	new.u<-old.u+rnorm(n)*usd
	new.u[1]<-0
	new.uvec<-rep(new.u,each=Tv)
	new.muvec<-X %*% old.beta + new.uvec
	new.res<-Y-new.muvec
	new.ssq<-(new.res)^2
	new.ilike<-apply(matrix(-0.5*new.ssq/old.sigma.e^2,nrow=n,byrow=T),1,sum)
	new.uilike<-dnorm(new.u,0,old.sigma.u,log=T)
	iacc<-log(runif(n)) < (new.ilike+new.uilike)-(old.ilike+old.uilike)
	old.u[iacc]<-new.u[iacc]
	old.uvec<-rep(old.u,each=Tv)	
	old.muvec<-new.muvec
	old.clike<--0.5*n*Tv*log(old.sigma.e^2)+sum(-0.5*(Y-old.muvec)^2)/old.sigma.e^2
	old.uilike[iacc]<-new.uilike[iacc]
	old.ulike<-sum(old.uilike[-1])

	#Update beta
	fb<-lm(Y~X-1+offset(old.uvec))
	beta.prop.m<-as.numeric(coef(fb))
	beta.prop.V<-summary(fb)$cov.unscaled
	new.beta<-as.numeric(mvrnorm(1,mu=beta.prop.m,Sigma=beta.prop.V))
	new.muvec<-X %*% new.beta + old.uvec
	new.clike<--0.5*n*Tv*log(old.sigma.e^2)+sum(-0.5*(Y-new.muvec)^2)/old.sigma.e^2
	new.prior.beta<--0.5*p*log(old.sigma.e^2*sigma.0^2)-0.5*sum(new.beta)^2/(old.sigma.e^2*sigma.0^2)
	old.term<-dmnorm(old.beta, mean = beta.prop.m, varcov=beta.prop.V, log = TRUE)
	new.term<-dmnorm(new.beta, mean = beta.prop.m, varcov=beta.prop.V, log = TRUE)

	if(log(runif(1)) < (new.clike+new.prior.beta-new.term)-(old.clike+old.prior.beta-old.term)){
		old.beta<-new.beta
		old.clike<-new.clike
		old.muvec<-X %*% old.beta + old.uvec
		old.res<-Y-old.muvec
		old.ssq<-(old.res)^2
		old.ilike<-apply(matrix(-0.5*old.ssq/old.sigma.e^2,nrow=n,byrow=T),1,sum)
	}

	#Update old.sigma.e
	old.sigma.e<-sqrt(1/rgamma(1,n*Tv/2+a.e,sum(old.ssq)/2+b.e))		
	old.clike<--0.5*n*Tv*log(old.sigma.e^2)-0.5*sum(old.ssq)/old.sigma.e^2
	old.prior.beta<--0.5*p*log(old.sigma.e^2*sigma.0^2)-0.5*sum(old.beta)^2/(old.sigma.e^2*sigma.0^2)

	#Update old.sigma.u
	new.sigma.u<-abs(rnorm(1,old.sigma.u,usd))
	new.uilike<-dnorm(old.u,0,new.sigma.u,log=T)
	new.ulike<-sum(new.uilike[-1])
	new.prior.sigu<-log(dinvgamma(new.sigma.u^2,a.u,b.u))

	if(log(runif(1))< (new.ulike+new.prior.sigu)-(old.ulike+old.prior.sigu)){
		old.sigma.u<-new.sigma.u
		old.uilike<-new.uilike
		old.ulike<-new.ulike
		old.prior.sigu<-new.prior.sigu
	}

	
	if(iter > nburn & iter %% nthin == 0){
		ico<-ico+1
		growth.beta[ico,]<-old.beta
		growth.sigmas[ico,]<-c(old.sigma.e,old.sigma.u)
		growth.us[ico,]<-old.u
		if(iter %% 50 == 0){print(c(iter,old.sigma.e,old.sigma.u))}

	}


}
(t1<-proc.time()[3]-t0)


bname<-c(expression(beta[1]),expression(beta[2]),expression(beta[3]))
pdf("Growth-Beta-Trace.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,2));for(j in 1:p){plot(growth.beta[,j],ylab=bname[j])}
dev.off()
sname<-c(expression(sigma[epsilon]),expression(sigma[upsilon]))
pdf("Growth-Sigma-Trace.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,1));for(j in 1:2){plot(growth.sigmas[,j],ylab=sname[j])}
dev.off()


uname<-c(expression(upsilon[2]),expression(upsilon[3]),expression(upsilon[4]),expression(upsilon[5]),expression(upsilon[6]))
pdf("Growth-Ups-Trace.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,3));for(j in 2:6){plot(growth.us[,j],ylab=uname[j],ylim=range(-0.5,0.5))}
dev.off()
pdf("Growth-Ups-Scatter.pdf",paper="USr",height=9,width=9)
par(mfrow=c(1,1));pairs(cbind(growth.us[,2:6],growth.sigmas[,2]),label=c(uname,expression(sigma[upsilon])))
dev.off()

ubname<-c(expression(upsilon[2]+beta[1]),expression(upsilon[3]+beta[1]),
expression(upsilon[4]+beta[1]),expression(upsilon[5]+beta[1]),expression(upsilon[6]+beta[1]))

pdf("Growth-Ups-ACF.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,3));for(j in 2:6){acf(growth.us[,j],main=uname[j])}
dev.off()

pdf("Growth-Ups+Beta-Trace.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,3));for(j in 2:6){plot(growth.us[,j]+growth.beta[,1],ylab=ubname[j],ylim=range(0.7,1.3))}
dev.off()

pdf("Growth-Ups+Beta-ACF.pdf",paper="USr",height=9,width=11)
par(mfrow=c(2,3));for(j in 2:6){acf(growth.us[,j]+growth.beta[,1],main=ubname[j])}
dev.off()

growth.ustars<-growth.us
for(j in 1:n){growth.ustars[,j]<-growth.ustars[,j]/growth.sigmas[,2]}

par(mfrow=c(1,1));pairs(growth.ustars[,2:6])

ub<-growth.us;for(j in 1:n){ub[,j]<-ub[,j]+growth.beta[,1]}

par(mfrow=c(1,1));pairs(ub[,2:6])