On Sun, 29 Nov 2015, Pindar Os wrote:
Hi there, I’m still working on the restricted OLS regressions and
in the meantime even uploaded an early version of a function package
that uses a SQP of GNU R and a bootstrap routine to get coeffs and
standard errors. In order to produce a nice show case of the
optimization problem I generate all combination of coefficients with
a step size of n. If there are 3 coeffs that are >=0 and should sum
up to 1 the matrix can be easily generated via loop. Since gretl has
a special "for loop" one could use this loop version in order to
fill the percentages in the matrix in one step. However, there is a
bug. Have a look. I’m using the current x64 snapshot for windows.
Cheers and a nice first Sunday in Advent Leon
<hansl>
scalar intN_1_a = 50 # step size of 2 percent
scalar intStep = 1/intN_1_a
eval intStep
matrix mCombis = NA
matrix mCombisALL = NA
set stopwatch
loop for i=0..intN_1_a --quiet
loop for j=0..intN_1_a --quiet
if i+j<=intN_1_a
if i+j == 0
mCombisALL = i ~ j ~ (intN_1_a-i-j)
else
mCombisALL |= i ~ j ~ (intN_1_a-i-j)
endif
endif
endloop
endloop
eval $stopwatch/60
eval rows(mCombisALL)
# correct number of rows since
scalar intSum3 = (intN_1_a)*3 + (1+intN_1_a-3+1)*(1+intN_1_a-3)/2
# for three coeefs I found this closed formula, with more it’s another…
eval intSum3
set stopwatch
loop for (i=0; i<=1; i+=intStep) --quiet
loop for (j=0; j<=1; j+=intStep) --quiet
if i+j<=1
if i+j == 0
mCombis = i ~ j ~ (1-i-j)
else
mCombis |= i ~ j ~ (1-i-j)
endif
endif
endloop
endloop
eval rows(mCombis)
eval $stopwatch
# correct number
scalar intSum3 = (intN_1_a)*3 + (1+intN_1_a-3+1)*(1+intN_1_a-3)/2
eval intSum3
# wrong, WHY? The reason has to do with the.xx999 values and the missing 1s.
<hansl>
It's discussed in the chapter on loops in the User's Guide. You need
to allow for the fact that when your i and j (in the second example)
"ought to" equal 1 if calculation were done with infinite precision,
they're likely to be slightly off. Hence you want something like
loop for (i=0; i<1.001; i+=intStep)
loop for (j=0; j<1.001; j+=intStep)
if i+j < 1.001
...
For reference, at the gretl console:
? printf "%.18g\n", 50*(1/50)
1.00000000000000022
Allin Cottrell