Count function calls

Thu, Mar 7, 2013 4:11 AM

Dear Giovanni,

apologize for this late reply! I was testing and reading a lot of stuff. I tried your suggestions and the problem of singularity in the regressor cross product vanishes when using the Group Mean function 'pgm' instead of 'pvcm'. Nevertheless, I found the collinearity in the regressors and 'pvcm' ran further, until I got another error. This time a little different but still directed towards the same area:

'Error in solve.default(crossprod(X[[i]])[!coefna[i, ], !coefna[i, ]])) :
	system is computationally singular: reciprocal condition number = 1.86131e-17'

Following the description of R's function 'rcond', the reciprocal condition number measures how close a matrix is to be rank deficient. So in this case one regressor crossproduct seems to be sufficiently close to singularity, such that inversion becomes impossible. 

I attached you a simple subsample. If you let run the following commands on it:

require(plm)
data <- read.csv(path, sep = ",", header = TRUE)
pdata <- plm.data(data, index = c("gvkey", "fyear"))
debug(plm:::pvcm)
model.pvcm <- pvcm(LOGDXSGA~LOGDSALE + DECRLOGDSALE + DECRLOGDSALEWDAVG, data = pdata, model = "random")
...... go until
Browse[2]> 
debug: ml <- split(data cond)

Then type:
A <- as.matrix(ml$'25062'[, 2:4])
XX <- t(A) %*% A
qr(XX)$rank
[1] 2                            Aha, the rank is only 2! 
rcond(XX) 
[1] 9.339817e-16    Alright, the matrix XX is pretty near to singularity, even it is not really singular!

Let us have a look at it:
A[, 2]/A[, 3]
.....                              A[, 2] is almost a multiple of A[, 3], if the latter is multiplied by -6.231979 or -6.231331 ..... at least it suffices to let the rank shrink to 2. 

The thing is, the same regression goes through in the Stata package. As I would say, I am an R fanatist, I would like to know, why it runs in Stata and not in R.....and that one can let it run in R as well. Different number formats/precision? Maybe it suffices in such a case to enforce full rank of the crossproduct by enforcing positive definitess, for instance via the function 'nearPD' in the R package 'Matrix':

Try 
library(Matrix)
qr(nearPD(XX)$mat)$rank
[1] 3                                                 Yey!..... :)

Tell me your opinion Giovanni! Give me a little help here please! 


Best 

Simon
 

P.S. Thank your for this fantastic package making panel data estimation possible in the R environment! 

-------------- next part --------------

On Feb 27, 2013, at 1:40 PM, Millo Giovanni <Giovanni_Millo at Generali.com> wrote:

Hello.
Another thing you may want to do depends on whether you are using
model="within" (the default) or model="random".

In the first case, pvcm() estimates separate regressions, so you just
need to loop lm() on individual indices to spot where it fails.

In the second case, what you may want to do is try a similar estimator:
pmg(..., model="mg"), which is an unweighted version of Swamy's
estimator in pvcm() and a simpler function to read and modify, possibly
using the global assignment operator '<<-' as already suggested by
William to output diagnostics.

Best,
Giovanni

-----Original Message-----
From: r-help-bounces at r-project.org

[mailto:r-help-bounces at r-project.org] On Behalf

Of Simon Zehnder
Sent: Tuesday, February 26, 2013 2:53 AM
To: r-help at r-project.org help
Subject: [R] Count function calls

Dear R-users,

I have the following problem: I am running the function 'pvcm' from

the 'plm' Panel Data

package. Inside this function 'solve' is called and gives for a

certain individual data series

an exception because of singularity. I would like to know which

individual data series

causes this error. I tried to debug it, but this is truly painful, as

solve is called inside of

'lapply' and there are over 5,000 individual data series in the panel.

Now, what I would like to do is to count the calls to 'solve' inside

the function, so I can

see, where the function throws the exception. I tried to use 'trace'

with a count variable,

but I have no clue how to define a global variable to be used by trace

and updated at

every call.....is there another approach?


Best

Simon

______________________________________________

Count function calls

Thread (2 messages)