[R] Quadratic programming

[Maija Sirkjärvi]

I really appreciate you helping me with this! I just don't seem to figure
it out.

(1) I don't know why you think bvec should be a matrix. The
documentation clearly says it should be a vector (implying not a
matrix).

- I've written it in a form of a matrix with one row and 2*J-3
columns. (0,2*J-3,1). I thought that it would pass as a vector. I made it a
vector now.

The thing is that I'm trying to replicate a C++ code with R. The C++ code
imposes shape restrictions on the function and works perfectly. The C++
code is attached and after that the same for R. You said you haven't used
the QP package for a decade. Is there a better/another package for these
types of problems?

Thanks again!

C++ code:

int main()
{
Print("Begin");

/* Bootstrap Parameters */
long RandomSeed = -98345; // Set random seed
const int S = 100; // Number of bootstrap replications

/* Housing Demand Parameters */
const double Beta =  1.1613;
const double v    =  0.7837;
const double Eta  = -0.5140;

/* Quadratic Programming Tolerance Parameters */
const double Delta1 =  0.000001;
const double Delta2 =  0.0001;

/* Read in Dataset */
DataFileDelim d("K:\\Data\\Local Jurisdictions\\AEJ Data.dat",'\t');
d.SortBy<double>("AfterTaxPrice");
int J = d.NumObs();
Vector<double> Price = d.Get<double>("AfterTaxPrice");
Vector<double> Educ  = d.Get<double>("EducAll");
Vector<double> Crime = d.Get<double>("CrimeIndex");
Vector<double> Dist  = d.Get<double>("RushTrav");

/* Adjust Crime Data */
for(int j=0; j<J; j++) if(Crime[j] < 0.005) Crime[j] = 0.0025;  // Set
crime to half of minimum community to avoid log(0)
for(int j=0; j<J; j++) Crime[j] = log(Crime[j]);

/* Compute Ranks */
Vector<double> Rank1(J);
for(int j=0; j<J; j++) Rank1[j] = (One + j) / (double)J;

/* Data for Non-parametric Regression */
Vector<double> X(J);
Matrix<double> Y(J,2);
Vector<double> Z(J);
for(int j=0; j<J; j++)
{
X(j)   = Rank1(j);
Y(j,0) = Crime(j);
Y(j,1) = Dist(j);
Z(j)   = Educ(j);
}

/* Constrained Estimates */
int RhoK1 = 12;
Vector<double> RhoGrid1 = Grid(-1.2,-0.1,RhoK1);
Matrix<double> gTrans(J,RhoK1);
for(int rk=0; rk<RhoK1; rk++)
{
/* Loop Over Rho */
double Rho = RhoGrid1(rk);

/* Print Progress */
cout << Rho << " ";
cout.flush();

/* Set Up Quadratic Programing Problem */
Vector<double> hSmooth(J);
for(int j=0; j<J; j++) hSmooth(j) = -pow(kr.Phi(j),Rho);
Vector<double> Q(J);
for(int j=0; j<J; j++) Q(j) = exp(-Rho * (Beta * pow(Price(j),Eta + One) -
One) / (One + Eta));
SymmetricMatrix<double> H(J,Zero);
Vector<double> c(J,Zero);
Matrix<double> Aeq(0,J);
Vector<double> beq(0);
Matrix<double> Aneq(2*J-3,J,Zero);
Vector<double> bneq(2*J-3);
Vector<double> lb(J,-Inf);
Vector<double> ub(J,Inf);
for(int j=0; j<J; j++) H(j,j) = One;
for(int j=0; j<J; j++) c(j) = -hSmooth(j);
for(int j=1; j<J; j++)
{
Aneq(j-1,j-1) = -One;
Aneq(j-1,j)   = One;
bneq[j-1]     = Delta1;
}
for(int j=2; j<J; j++)
{
Aneq(J-1+j-2,j)   = -One / (Q(j) - Q(j-1));
Aneq(J-1+j-2,j-1) = One / (Q(j) - Q(j-1)) + One / (Q(j-1) - Q(j-2));
Aneq(J-1+j-2,j-2) = -One / (Q(j-1) - Q(j-2));
bneq[J-1+j-2]     = Delta2;
}

/* Solve Constrained Optimization Problem Using Quadratic Programming */
MinQuadProg qp(c,H,Aeq,beq,Aneq,bneq,lb,ub);
qp.PrintLevel = 0;
qp.Solve();

/* Constrained Estimate */
for(int j=0; j<J; j++) gTrans(j,rk) = pow(-qp.Solution(j),One / Rho);


And my version for R:

Beta =  1.1613;
v    =  0.7837;
Eta  = -0.5140;
Delta1 =  0.000001;
Delta2 =  0.0001;

newdata <- AEJData[order(AEJData$AfterTaxPrice),]
View(newdata)
Price = (newdata$AfterTaxPrice)
Educ  = (newdata$EducAll)
Crime = (newdata$CrimeIndex)
Dist  = (newdata$RushTrav)

install.packages("sp")
install.packages("rgdal")
install.packages("raster")
install.packages("calibrate")
install.packages("zeros")
install.packages("quadprog")
library(sp)
library(rgdal)
library(raster)
library(calibrate)
library(quadprog)


J <- length(Price)
hs <- numeric(J)
for(j in 1:J){
  hs[j] <-(-(gEst$KernelRegPartLin..Phi[j]^(-0.1)))
}
hs

Q <- rep(0,J)
for(j in 1:(length(Price))){
  Q[j] <- exp((-0.1) * (Beta *Price[j]^(Eta + 1) - 1) / (1 + Eta))
}
Q
plot(Q)


Dmat <- matrix(0,nrow= J, ncol=J)
diag(Dmat) <- 1
dvec <- -hs
Aeq <- 0
beq <- 0
Amat <- matrix(0,J,2*J-3)
bvec <- rep(0,2*J-3)

bv <- t(bvec)
bv

for(j in 2:nrow(Amat)){
  Amat[j-1,j-1] = -1
  Amat[j,j-1] = 1
}

for(j in 3:nrow(Amat)){
  Amat[j,J+j-3] = -1/(Q[j]-Q[j-1])
  Amat[j-1,J+j-3] = 1/(Q[j]-Q[j-1])
  Amat[j-2,J+j-3] = -1/(Q[j-1]-Q[j-2])
}

for(j in 2:ncol(bv)) {
  bv[j-1] = Delta1
}
for(j in 3:ncol(bv)) {
  bv[J-1+j-2] = Delta2
}

solution <- solve.QP(Dmat,dvec,Amat,bvec)

solution

gt <- matrix(0,J,1)
gt

Rho <- c(-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1)

Rh <- t(Rho)

gt= (-solution)^(1/(Rh))
gt


ma 21. syysk. 2020 klo 23.38 Abby Spurdle (spurdle.a using gmail.com) kirjoitti:

> I was wondering if you're trying to fit a curve, subject to
> monotonicity/convexity constraints...
> If you are, this is a challenging topic, best of luck...
>
>
> On Tue, Sep 22, 2020 at 8:12 AM Abby Spurdle <spurdle.a using gmail.com> wrote:
> >
> > Hi,
> >
> > Sorry, for my rushed responses, last night.
> > (Shouldn't post when I'm about to log out).
> >
> > I haven't used the quadprog package for nearly a decade.
> > And I was hoping that an expert using optimization in finance in
> > economics would reply.
> >
> > Some comments:
> > (1) I don't know why you think bvec should be a matrix. The
> > documentation clearly says it should be a vector (implying not a
> > matrix).
> > The only arguments that should be matrices are Dmat and Amat.
> > (2) I'm having some difficulty following your quadratic program, even
> > after rendering it.
> > Perhaps you could rewrite your expressions, in a form that is
> > consistent with the input to solve.QP. That's a math problem, not an R
> > programming problem, as such.
> > (3) If that fails, then you'll need to produce a minimal reproducible
> example.
> > I strongly recommend that the R code matches the quadratic program, as
> > closely as possible.
> >
> >
> > On Mon, Sep 21, 2020 at 9:28 PM Maija Sirkjärvi
> > <maija.sirkjarvi using gmail.com> wrote:
> > >
> > > Hi!
> > >
> > > I was wondering if someone could help me out. I'm minimizing a
> following
> > > function:
> > >
> > > \begin{equation}
> > > $$\sum_{j=1}^{J}(m_{j} -\hat{m_{j}})^2,$$
> > > \text{subject to}
> > > $$m_{j-1}\leq m_{j}-\delta_{1}$$
> > > $$\frac{1}{Q_{j-1}-Q_{j-2}} (m_{j-2}-m_{j-1}) \leq
> \frac{1}{Q_{j}-Q_{j-1}}
> > > (m_{j-1}-m_{j})-\delta_{2} $$
> > > \end{equation}
> > >
> > > I have tried quadratic programming, but something is off. Does anyone
> have
> > > an idea how to approach this?
> > >
> > > Thanks in advance!
> > >
> > > Q <- rep(0,J)
> > > for(j in 1:(length(Price))){
> > >   Q[j] <- exp((-0.1) * (Beta *Price[j]^(Eta + 1) - 1) / (1 + Eta))
> > > }
> > >
> > > Dmat <- matrix(0,nrow= J, ncol=J)
> > > diag(Dmat) <- 1
> > > dvec <- -hs
> > > Aeq <- 0
> > > beq <- 0
> > > Amat <- matrix(0,J,2*J-3)
> > > bvec <- matrix(0,2*J-3,1)
> > >
> > > for(j in 2:nrow(Amat)){
> > >   Amat[j-1,j-1] = -1
> > >   Amat[j,j-1] = 1
> > > }
> > > for(j in 3:nrow(Amat)){
> > >   Amat[j,J+j-3] = -1/(Q[j]-Q[j-1])
> > >   Amat[j-1,J+j-3] = 1/(Q[j]-Q[j-1])
> > >   Amat[j-2,J+j-3] = -1/(Q[j-1]-Q[j-2])
> > > }
> > > for(j in 2:ncol(bvec)) {
> > >   bvec[j-1] = Delta1
> > > }
> > > for(j in 3:ncol(bvec)) {
> > >   bvec[J-1+j-2] = Delta2
> > > }
> > > solution <- solve.QP(Dmat,dvec,Amat,bvec=bvec)
> > >
> > >         [[alternative HTML version deleted]]
> > >
> > > ______________________________________________
> > > R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
>

	[[alternative HTML version deleted]]