[R] Restricted Least Squares

Simon Wood snw at mcs.st-and.ac.uk
Tue Apr 9 12:10:19 CEST 2002


> I need help regarding estimating a linear model where restrictions are imposed on the coefficients. An example is as follows:
> Y_{t+2}=a1Y_{t+1} + a2 Y_t + b x_t + e_t
> restriction
> a1+ a2 =1
> Is there a function or a package that can estimate the coefficient of a model like this? I want to estimate the coefficients rather than test them.
- pcls() in package mgcv will do this (although you'll need to
re-parameterize so that the restriction is a1+a2=0). 

Alternatively the following fits: 
y = b0 + b1 x + e subject to b0 + b1 = 1
having reparameterized as
y = b0 + (b1+1) x + e subject to b0 + b1 = 0

x<-runif(100)
y<-x+rnorm(100)*0.1            # simulate some data
cm<-matrix(c(1,1),1,2)         # make constraint matrix C so Cb=0 
cqr<-qr(t(cm))                 # QR decomp of C'
# if m constraints then last m cols of Q give null space Z of constraints: 
Z<-qr.Q(cqr,complete=TRUE)[,2] 
X<-model.matrix(y~x)            # model matrix
XZ<-X%*%Z                     # model matrix in constraint null space
p<-lm(y~XZ-1+offset(x))$coeff # estimate model in null space 
p<-Z%*%as.matrix(p)           # estimates in original space 
p[2]<-p[2]+1                  # undo reparameterization
p      
plot(x,y);abline(p)

... this can be easily generalized to do what you want.

Simon

  ______________________________________________________________________
> Simon Wood  snw at st-and.ac.uk  http://www.ruwpa.st-and.ac.uk/simon.html
> The Mathematical Institute, North Haugh, St. Andrews, Fife KY16 9SS UK
> Direct telephone: (0)1334 463799          Indirect fax: (0)1334 463748 




-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._



More information about the R-help mailing list