Quadratic programming

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

Quadratic programming

Maija Sirkjärvi
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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Quadratic programming

aBBy Spurdle, ⍺XY
Are you using the quadprog package?
If I can take a random shot in the dark, should bvec be -bvec?


On Mon, Sep 21, 2020 at 9:28 PM Maija Sirkjärvi
<[hidden email]> 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]]
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Quadratic programming

aBBy Spurdle, ⍺XY
Sorry, ignore the last part.
What I should have said, is the inequality has the opposite sign.
>= bvec (not <= bvec)


On Mon, Sep 21, 2020 at 10:05 PM Abby Spurdle <[hidden email]> wrote:

>
> Are you using the quadprog package?
> If I can take a random shot in the dark, should bvec be -bvec?
>
>
> On Mon, Sep 21, 2020 at 9:28 PM Maija Sirkjärvi
> <[hidden email]> 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]]
> >
> > ______________________________________________
> > [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Quadratic programming

aBBy Spurdle, ⍺XY
One more thing, is bvec supposed to be a matrix?

Note you may need to provide a reproducible example, for better help...

On Mon, Sep 21, 2020 at 10:09 PM Abby Spurdle <[hidden email]> wrote:

>
> Sorry, ignore the last part.
> What I should have said, is the inequality has the opposite sign.
> >= bvec (not <= bvec)
>
>
> On Mon, Sep 21, 2020 at 10:05 PM Abby Spurdle <[hidden email]> wrote:
> >
> > Are you using the quadprog package?
> > If I can take a random shot in the dark, should bvec be -bvec?
> >
> >
> > On Mon, Sep 21, 2020 at 9:28 PM Maija Sirkjärvi
> > <[hidden email]> 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]]
> > >
> > > ______________________________________________
> > > [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Quadratic programming

Maija Sirkjärvi
Thank you for your response!

Bvec is supposed to be a matxit. I'm following the solve.QP (
https://www.rdocumentation.org/packages/quadprog/versions/1.5-8/topics/solve.QP
).

I'm not sure what would be the best way to solve a quadratic programming
problem in R.


ma 21. syysk. 2020 klo 13.20 Abby Spurdle ([hidden email]) kirjoitti:

> One more thing, is bvec supposed to be a matrix?
>
> Note you may need to provide a reproducible example, for better help...
>
> On Mon, Sep 21, 2020 at 10:09 PM Abby Spurdle <[hidden email]> wrote:
> >
> > Sorry, ignore the last part.
> > What I should have said, is the inequality has the opposite sign.
> > >= bvec (not <= bvec)
> >
> >
> > On Mon, Sep 21, 2020 at 10:05 PM Abby Spurdle <[hidden email]>
> wrote:
> > >
> > > Are you using the quadprog package?
> > > If I can take a random shot in the dark, should bvec be -bvec?
> > >
> > >
> > > On Mon, Sep 21, 2020 at 9:28 PM Maija Sirkjärvi
> > > <[hidden email]> 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]]
> > > >
> > > > ______________________________________________
> > > > [hidden email] 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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Quadratic programming

aBBy Spurdle, ⍺XY
In reply to this post by Maija Sirkjärvi
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
<[hidden email]> 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]]
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Quadratic programming

aBBy Spurdle, ⍺XY
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 <[hidden email]> 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
> <[hidden email]> 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]]
> >
> > ______________________________________________
> > [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: 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 ([hidden email]) 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 <[hidden email]> 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
> > <[hidden email]> 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]]
> > >
> > > ______________________________________________
> > > [hidden email] 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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Quadratic programming

aBBy Spurdle, ⍺XY
> I'm trying to replicate a C++ code with R.

Notes:
(1) I'd recommend you make the code more modular.
i.e. One function for initial data prep/modelling, one function for
setting up and solving the QP, etc.
This should be easier to debug.
(However, you would probably have to do it to the C++ code first).
(2) Your R code is not completely reproducible.
i.e. AEJData
(3) For the purposes of a reproducible example, your code can be simplified.
i.e. Only one contributed R package should be attached.

Regardless of (1) above, you should be able to identify at what point
the C++ and R code becomes inconsistent.
The simplest approach is to add print-based functions into both the
C++ and R code, and print out state data, at each major step.
Then all you need to do is compare the output for both.

> Is there a better/another package for these types of problems?

I'm not sure.
After a quick search, this is the best I found:

scam::scam
scam::shape.constrained.smooth.terms

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Quadratic programming

Maija Sirkjärvi
Thank you for giving me your time!

The problem is the quadratic optimization part. Something goes wrong along
the way. In C++ loops run from 0 and in R they run from 1, and I've tried
to take that into account. Still I'm having hard time figuring out the
mistake I make, cause I get a result from my R code. It's just not the same
that I get with the C++.

Here are the quadratic optimization parts for both codes.

C++

  /* 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();

And R

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)

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:nrow(bvec)) {
  bvec[j-1] = Delta1
}
for(j in 3:nrow(bvec)) {
  bvec[J-1+j-2] = Delta2
}

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





ke 23. syysk. 2020 klo 9.12 Abby Spurdle ([hidden email]) kirjoitti:

> > I'm trying to replicate a C++ code with R.
>
> Notes:
> (1) I'd recommend you make the code more modular.
> i.e. One function for initial data prep/modelling, one function for
> setting up and solving the QP, etc.
> This should be easier to debug.
> (However, you would probably have to do it to the C++ code first).
> (2) Your R code is not completely reproducible.
> i.e. AEJData
> (3) For the purposes of a reproducible example, your code can be
> simplified.
> i.e. Only one contributed R package should be attached.
>
> Regardless of (1) above, you should be able to identify at what point
> the C++ and R code becomes inconsistent.
> The simplest approach is to add print-based functions into both the
> C++ and R code, and print out state data, at each major step.
> Then all you need to do is compare the output for both.
>
> > Is there a better/another package for these types of problems?
>
> I'm not sure.
> After a quick search, this is the best I found:
>
> scam::scam
> scam::shape.constrained.smooth.terms
>

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Quadratic programming

aBBy Spurdle, ⍺XY
Before going any further, I have to check, what is:
MinQuadProg qp

Also, if I'm following the C++ code correctly, H, is an identity matrix.
This implies the input to the C++ solver, requires the QP in a
different form to the R solver.
In which case, the C++ inputs and the R inputs, should be different...?

(A)
It may be worthwhile comparing the solver output (for C++ and R) using
a *much* simpler example.
e.g. the examples from the quadprog package.

(B)
If they produce the same output (using a simple example), then that
implies there's a difference in your inputs.
So, you just need to work out which input values are different.
Expanding on my previous post, just print them out.
But check (A) above first.









On Thu, Sep 24, 2020 at 11:15 PM Maija Sirkjärvi
<[hidden email]> wrote:

>
> Thank you for giving me your time!
>
> The problem is the quadratic optimization part. Something goes wrong along the way. In C++ loops run from 0 and in R they run from 1, and I've tried to take that into account. Still I'm having hard time figuring out the mistake I make, cause I get a result from my R code. It's just not the same that I get with the C++.
>
> Here are the quadratic optimization parts for both codes.
>
> C++
>
>   /* 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();
>
> And R
>
> 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)
>
> 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:nrow(bvec)) {
>   bvec[j-1] = Delta1
> }
> for(j in 3:nrow(bvec)) {
>   bvec[J-1+j-2] = Delta2
> }
>
> solution <- solve.QP(Dmat,dvec,Amat,bvec)
>
>
>
>
>
> ke 23. syysk. 2020 klo 9.12 Abby Spurdle ([hidden email]) kirjoitti:
>>
>> > I'm trying to replicate a C++ code with R.
>>
>> Notes:
>> (1) I'd recommend you make the code more modular.
>> i.e. One function for initial data prep/modelling, one function for
>> setting up and solving the QP, etc.
>> This should be easier to debug.
>> (However, you would probably have to do it to the C++ code first).
>> (2) Your R code is not completely reproducible.
>> i.e. AEJData
>> (3) For the purposes of a reproducible example, your code can be simplified.
>> i.e. Only one contributed R package should be attached.
>>
>> Regardless of (1) above, you should be able to identify at what point
>> the C++ and R code becomes inconsistent.
>> The simplest approach is to add print-based functions into both the
>> C++ and R code, and print out state data, at each major step.
>> Then all you need to do is compare the output for both.
>>
>> > Is there a better/another package for these types of problems?
>>
>> I'm not sure.
>> After a quick search, this is the best I found:
>>
>> scam::scam
>> scam::shape.constrained.smooth.terms

______________________________________________
[hidden email] 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.