R help - Sep 2008 - panel data analysis possible with mle2 (bbmle)?

Dear R community,

I want to estimate coefficients in a (non-linear) system of equations using
'mle2' from the "bbmle" package. Right now the whole data is
read in as just
one long time series, when it's actually 9 cross sections with 30
observations
each. I would like to be able to test and correct for autocorrelation but
haven't found a way to do this in this package. 
On the other hand I have found a package called "plm" serving for
panel
likelihood maximization, but here the regression formula apparently must be
one single equation, whereas I have a system of equations.

In order to give you an idea of what I'm trying to do, I have written this
-as
short as possible- sketch of the model - sorry I couldn't make it any
shorter:

################################################################
cross <- 3;years <- 13 #cross sections and years to be included

frontr<-vector(length=years);ggwth <- 0.03; frontr[1]=17
for (i in 2:years) {frontr[i]=frontr[i-1]*(1+ggwth)}
frontr <- rep(frontr,times=cross)

atech <- matrix(0,years,cross);rr1 <- matrix(0,years,cross)
agwth <- matrix(0,years,cross);cfgwth <- matrix(0,years,cross)

aobs <-
c(15.89,16.69,17.93,17.49,18.3,19.1,19.2,19.4,20.29,21.13,21.42,22.76,23.83,
         
14.1,14.4,13.4,14.9,15.23,15.4,15.55,16.7,17.8,18.87,18.99,19.24,20.59,
          14.3,14.4,14.7,14.9,15.2,15.43,15.5,17.1,17.5,18,18.3,18.7,19);
dim(aobs) <- c(years,cross)

ra <- c(-0.62, -0.81, -1.05, -0.97, -0.87, -0.66, -0.67, -0.56,
    -0.333, -0.219, -0.13, -0.116, -0.474, -0.838, -0.821, -0.95,
    -1.03, -1.26, -1.21, -1.12, -1.19, -1.14, -1, -0.97,
    -0.92, -1.05, 0.04, 0.02, -0.11, -0.36, -0.57, -0.7,
    -0.71, -0.61, -0.68, -0.9, -1, -1.04, -1.1)

sa <- c(1.79, 1.81, 1.81, 1.81, 1.81, 1.81, 1.81, 1.81, 1.81,
    1.82, 1.8, 1.79, 1.783, 0.93, 0.96, 0.99, 1.02, 1.04, 1.07,
    1.1, 1.12, 1.14, 1.17, 1.19, 1.21, 1.23, 1.35, 1.4,
    1.45, 1.5, 1.56, 1.61, 1.66, 1.71, 1.75, 1.81, 1.85,
    1.89, 1.93)

ta <- c(-1.06, -1.01, -0.97, -0.94, -0.89, -0.82, -0.75, -0.69,
    -0.58, -0.54, -0.49, -0.39, -0.31, 0.02, 0.07, 0.087,
    0.138, 0.132, 0.244, 0.354, 0.421, 0.606, 0.74, 0.818, 1.048,
    1.229, -1.17, -1.159, -1.16, -1.143, -1.137, -1.085, -1.028,
    -0.965, -0.919, -0.942, -0.902, -0.828, -0.79)
################this is the function that is called by mle2
exfunc <-
function(c1,c2,c3,cs1,cs2,cs3,ic1,ic2,ic3,cc1,cc2,cc3,alstar,beta1){
x<-matrix(0,years*cross,3)
   x[,1] <- 1
   x[,2] <- ra
   x[,3] <- sa
coeffs <- vector(length=3)
   coeffs[1] <- c1
   coeffs[2] <- c2
   coeffs[3] <- c3

csp <- rep(c(cs1,cs2,cs3),each=years)
e1 <- (x %*% coeffs);frfact <- csp*plogis(e1);cfrontr <- frfact*frontr;
cfrontr1 <- cfrontr;dim(cfrontr1) <- c(years,cross)
cfgwth[2:years,] <- (cfrontr1[2:years,]/cfrontr1[1:(years-1),])-1

x<-matrix(1,years*cross,3)
   x[,1] <- 1
   x[,2] <- sa
   x[,3] <- ta
coeffs <- vector(length=3)
   coeffs[1] <- cc1
   coeffs[2] <- cc2
   coeffs[3] <- cc3

e1 <- (x %*% coeffs);alpha1 <- alstar*plogis(e1);dim(alpha1) <-
c(years,cross)

icfr <- c(ic1,ic2,ic3)
atech[1,1:cross] <- icfr[1:cross]*cfrontr1[1,1:cross]
rr1[1,1:cross] <- atech[1,1:cross]/cfrontr1[1,1:cross]

for (i in 2:years) {
     
agwth[i,1:cross]=(rr1[i-1,1:cross]-alpha1[i-1,1:cross])*beta1*(1-rr1[i-1,1:cross])+cfgwth[i,1:cross]
      atech[i,1:cross]=atech[i-1,1:cross]*(1+agwth[i,1:cross])
      rr1[i,1:cross]=atech[i,1:cross]/cfrontr1[i,1:cross]}

pcdiffun <- function(a,b) (a-b)/b
pcdif <- pcdiffun(atech[2:years,],aobs[2:years,]);sd1 = sd(pcdif)
rt1 = -sum(dnorm(pcdif,mean=0,sd=sd1,log=TRUE))}
##############here is where the function ends 
library(bbmle)
maxlk <- mle2(exfunc,
              start=list(c1 = -3.74466306894815,c2 = -0.240135942405391,c3
4.44168922065451,
              cs1 = 0.962506678868304,cs2 = 0.965725991431272,cs3
0.80188688385734,
              ic1 = 2.25406162804465,ic2 = 1.68235691305909,ic3
2.13468107737406,
              cc1 = -8.564358,cc2 = 10.7195144315809,cc3 = 0,
              alstar = -2.59819924393425,beta1 = 0.0937403319000222),
              method="BFGS",control=list(maxit=1000))
summary(maxlk)
################################################################

The question now is, if anybody knows a way to make a panel data analysis with
'mle2' or alternatively, how to enter a system of equations instead of
just a
single formula in 'plm'.

Thanks in advance,
Erich