Or there is no way to do that in MATLAB..? the stats toolbox is pretty weakAlso, how to do regression with the autoregressive residuals?y(t) = b*x(t) + e(t)e(t) = a*e(t-1)+...plz advise, thanks!

MATLAB is generally pretty weak when it comes to providing canned routines. If you want software that would solve your problem with a single command, look into something like gretl, STATA, or EVIEWS. As you know, the formula for GLS is trivial: B = inv(X'inv(G)X)X'inv(G)yPerhaps Matlab has it built it, but otherwise you can use something like B =solve(inv(X'*inv(G)*X)*X'*inv(G),y) (typing the first formula directly is bad for numerical accuracy).Of course, you need to know what goes inside of G, which is the covariance matrix of e_0,e_1,...,e_T. For AR(1) model it has a trivial structure where all off diagonal entries depend on a single parameter that dies out geometrically. You get the parameter by estimating the AR model first, then building G, and solving the GLS equations above. Any good econometrics text book works through this example. I suggest that you look into Goldberger's graduate econometrics book, the chapter on GLS with serial correlation. Alternatively, forget the GLS, and do OLS with serial correlation robust standard errors. (again, MATLAB probably does not have functions that compute robust standard errors)

MATLAB does have it built in look up "rdivide"the logic is that if y=x*b +errorb(minimum (least square) error) b= x\y(y, x are in columns , the first row is the first observation, t=1, so the dimensions are y(MxN) , x(MxK) , and B(KxN))

You're right about the X\y trick, but that just solves the OLS problem. It would also solve the GLS problem after the X and y matrices are pre-multiplied by a Cholesky or some other factor of inv(G). It turns out there is no point in doing that as MATLAB has the lscov function for GLS in the presence of the known covariance matrix. But one still needs to find the covariance matrix for the AR(1) errors.One way of doing this is to compute OLS equation first, extract the residuals, and then fit AR(1) model on them. Then it is easy to build the G matrix because for the AR(1) model: e_t = B * e_t-1 + w_t,cov[e_t,e_(t-j)] = Var(e_t) * B^jThe constant Var(e_t) factor can be dropped from G because it will cancel out through the GLS equation.Any econometrics text has the implementation details for this procedure.