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 (spurdle.a at gmail.com) 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