# powers to pass to fdjac  p(x) is x   p(x,3) is x^3

function scalar p(scalar x, scalar p[1])

    return x^p

end function

function matrix fun4(strings fun, matrix rhs, matrix beta, matrix vcov, scalar s, scalar df)

    b = beta

    coll = rows(b)

    r = nelem(fun)

    matrix matr = zeros(r,coll)

    matrix matr2 = zeros(r,1)

    loop i=1..r --quiet

        fu = fun[i]

        vect = fdjac(b,@fu)

        matr[i,] = vect

        matr2[i,] = @fu-rhs[i]

    endloop

    ret = matr*vcov*transp(matr)

    reti = inv(ret)

    ff = transp(matr2)*reti*matr2/r

    pval = 1-cdf(F,r,df,ff)

    out = {ff,pval}

    colnames(out, "F p-value")

    return out

end function

## examples of restrictions

matrix rhs = {0,0,0}

strings S = array(2)

S[1] = "p(b[2])+p(b[3])"

S[2] = "p(b[2])+2*p(b[3])"

strings S1 = array(1)

S1[1] = "p(b[2])+p(b[3])"

open denmark.gdt

ols LRM const LRY IBO

eval fun4(S1,rhs,$coeff,$vcv,$sigma,$df)

eval fun4(S,rhs,$coeff,$vcv,$sigma.$df)