library(sn)
library(mvtnorm)
library(normalp)

set.seed(20120301)

### normal data
mL<-5
mS<-2
vL<-2
vS<-1
vLS<-0
vE<-1

TT<-5
N<-500
R<-100

mu<-c(mL,mS)
sigma<-array(c(vL,vLS,vLS,vS), dim=c(2,2))


for (k in 1:R){
y<-array(NA, dim=c(N, TT))
LS<-array(NA, dim=c(N, 2))

error<-NULL

for (t in 1:TT){
	error<-cbind(error, rnorm(N,df=df)*sqrt(vE))
}

for (i in 1:N){
  LS[i,]<-rmvnorm(1, mu, sigma)
  for (j in 1:TT){
    y[i, j] <- LS[i,1] + LS[i,2]*j + error[i,j]
  }
}


filename<-paste('data/data',k, '.txt',sep='')
write.table(y,filename,row.names=F,col.names=F)
}

### t data
mL<-5
mS<-2
vL<-2
vS<-1
vLS<-0
vE<-1/3
df<-3

TT<-5
N<-500
R<-100

mu<-c(mL,mS)
sigma<-array(c(vL,vLS,vLS,vS), dim=c(2,2))


for (k in 1:R){
y<-array(NA, dim=c(N, TT))
LS<-array(NA, dim=c(N, 2))

error<-NULL

for (t in 1:TT){
	error<-cbind(error, rt(N,df=df)*sqrt(vE))
}

for (i in 1:N){
  LS[i,]<-rmvnorm(1, mu, sigma)
  for (j in 1:TT){
    y[i, j] <- LS[i,1] + LS[i,2]*j + error[i,j]
  }
}


filename<-paste('data/data',k, '.txt',sep='')
write.table(y,filename,row.names=F,col.names=F)
}

### Exponential power data
sig<-c(1,1,1,1)
beta<- c(-.8, -.8, .8, .8) #rep(-.8, 5)
mL<-5
mS<-2
vL<-2
vS<-1
vLS<-0

TT<-4
N<-500

p<-2/(1+beta)
s<-sig*sqrt(gamma(1/p)/gamma(3/p))/p^(1/p)


mu<-c(mL,mS)
sigma<-array(c(vL,vLS,vLS,vS), dim=c(2,2))

y<-array(NA, dim=c(N, TT))
LS<-array(NA, dim=c(N, 2))

for (k in 1:100){
error<-NULL

for (t in 1:TT){
	error<-cbind(error, rnormp(N,0,s[t],p[t]))
}

for (i in 1:N){
  LS[i,]<-rmvnorm(1, mu, sigma)
  for (j in 1:TT){
    y[i, j] <- LS[i,1] + LS[i,2]*j + error[i,j]
  }
}

filename<-paste('data/t4data',k, '.txt',sep='')
write.table(y,filename,row.names=F,col.names=F)
}

## T=5
sig<-c(1,1,1,1,1)
beta<- c(-.8, -.8, 0, .8, .8) #rep(-.8, 5)
mL<-5
mS<-2
vL<-2
vS<-1
vLS<-0

TT<-5
N<-500

p<-2/(1+beta)
s<-sig*sqrt(gamma(1/p)/gamma(3/p))/p^(1/p)


mu<-c(mL,mS)
sigma<-array(c(vL,vLS,vLS,vS), dim=c(2,2))

y<-array(NA, dim=c(N, TT))
LS<-array(NA, dim=c(N, 2))

for (k in 1:100){
error<-NULL

for (t in 1:TT){
	error<-cbind(error, rnormp(N,0,s[t],p[t]))
}

for (i in 1:N){
  LS[i,]<-rmvnorm(1, mu, sigma)
  for (j in 1:TT){
    y[i, j] <- LS[i,1] + LS[i,2]*j + error[i,j]
  }
}

filename<-paste('data/t5data',k, '.txt',sep='')
write.table(y,filename,row.names=F,col.names=F)
}


### Skew normal data
mL<-5
mS<-2
vL<-2
vS<-1
vLS<-0
vE<-sqrt(pi/(pi-1))
alpha <- 10
delta<-alpha/sqrt(1+alpha^2)
omega<-sqrt(1/(1-2*delta^2/pi))
s.mu<- -omega*delta*sqrt(2/pi)

TT<-5
N<-500

mu<-c(mL,mS)
sigma<-array(c(vL,vLS,vLS,vS), dim=c(2,2))

y<-array(NA, dim=c(N, TT))
LS<-array(NA, dim=c(N, 2))
my<-array(NA, dim=c(N, TT))

for (k in 1:100){
error<-NULL

for (t in 1:TT){
	error<-cbind(error, rsn(N,s.mu, omega, alpha))
}

for (i in 1:N){
  LS[i,]<-rmvnorm(1, mu, sigma)
  for (j in 1:TT){
  	my[i,j]<- LS[i,1] + LS[i,2]*j
    #y[i, j] <- rsn(1,my[i,j],vE,alpha) #+ error[i,j]
	y[i, j] <- my[i,j] + error[i,j]
  }
}

filename<-paste('data/data',k, '.txt',sep='')
write.table(y,filename,row.names=F,col.names=F)
}