Re: [R] Calling stargazer() with do.call() in R 4.2.0
Dear Andrew
Thanks a lot for investigating the problem and for suggesting a
solution and a workaround! I have asked the maintainer of the
'stargazer' package to fix the problem. Until then, we can use the
workaround that you suggested.
Best regards,
Arne
On Sat, 28 May 2022 at 23:13, Andrew Simmons wrote:
>
> Hello,
>
>
> I don't have the slightest clue what stargazer is supposed to be
> doing, but it seems as though it's trying to create names for the ...
> list using:
> object.names.string <- deparse(substitute(list(...)))
>
> It makes an error in assuming that the return value of deparse will be
> a character string. For stargazer(res), object.names.string becomes:
> "list(res)"
>
> while for do.call(stargazer, list(res)), object.names.string becomes:
> [1] "list(structure(list(coefficients = c(`(Intercept)` =
> 6.41594246095523, "
> [2] "UrbanPop = 0.0209346588197249), residuals = c(Alabama =
> 5.56984732750071, "
> [3] "Alaska = 2.57919391569797, Arizona = 0.009284833466778, Arkansas
> = 1.33732459805853, "
> [4] "California = 0.679003586449805, Colorado = -0.148845848893771, "
> [5] "Connecticut = -4.72791119007405, Delaware = -2.02323789597542, "
> [6] "Florida = 7.30928483346678, Georgia = 9.72797800986128, Hawaii =
> -2.8535191429924, "
> [7] "Idaho = -4.94641403722037, Illinois = 2.2464808570076, Indiana =
> -0.576695284237348, "
> [8] "Iowa = -5.40921801367955, Kansas = -1.79762994305707, Kentucky =
> 2.19545528041907, "
> [9] "Louisiana = 7.60237005694293, Maine = -5.3836100607612, Maryland
> = 3.4814353981232, "
> [ reached getOption("max.print") -- omitted 110 entries ]
>
> perhaps the package maintainer could change the line from:
> object.names.string <- deparse(substitute(list(...)))
> to:
> object.names.string <- deparse1(substitute(list(...)), collapse = "")
>
> or you could change your code to:
> do.call(stargazer, alist(res))
>
> please note that using alist instead of list is only a workaround, you
> should still let the package maintainer know of this bug. If the
> maintainer asks, this is what I used to get the strings above:
> fun <- \(...) deparse(substitute(list(...)))
> data("USArrests")
> res <- lm( Murder ~ UrbanPop, data = USArrests)
> fun(res)
> print(do.call("fun", list(res)), max = 9)
>
> On Sat, May 28, 2022 at 4:41 PM Arne Henningsen
> wrote:
> >
> > On Sat, 28 May 2022 at 01:21, Uwe Ligges
> > wrote:
> > > On 27.05.2022 17:29, Arne Henningsen wrote:
> > >> Dear all (cc Marek = maintainer of the stargazer package)
> > >>
> > >> We use do.call() to automatically create many LaTeX tables with
> > >> stargazer but after upgrading to R 4.2.0, this no longer works. I
> > >> illustrate this with a simple reproducible example:
> > >>
> > >> R> data("USArrests")
> > >> R> res <- lm( Murder ~ UrbanPop, data = USArrests )
> > >> R> library(stargazer)
> > >> R> stargazer(res) # works as expected
> > >> R> do.call( stargazer, list(res) )
> > >> Error in if (is.na(s)) { : the condition has length > 1
> > >
> > > Without looking at the code in detail: The line aboce suggests the code
> > > needs an any():if(any(is.na(x))) raher than if(is.na(x)).
> >
> > Yes, this is likely a problem in the 'stargazer' package.
> >
> > ... but why does the problem occur when using do.call( stargazer, )
> > but the problem does *not* occur when using stargazer() directly?
> >
> > Best regards,
> > Arne
> >
> > >> Any ideas what we can do so that the last command works with R 4.2.0?
> > >>
> > >> /Arne
> >
> > --
> > Arne Henningsen
> > http://www.arne-henningsen.name
> >
> > __
> > R-help@r-project.org 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.
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org 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.
Re: [R] Interpreting coefficient in selection and outcome Heckman models in sampleSelection
Dear Marinella Implementing a function in the sampleSelection package that calculates marginal effects is still on my to-do list but I probably won't implement it soon, because I have many other things with higher priority. Sorry! However, you are invited to implement this feature in the sampleSelection package; I would assist you with this. Best regards, Arne On Mon, 28 Dec 2020 at 15:22, Marinella Cirillo via R-help wrote: > > Dear Arne, > > I have just read the exchange of messages with Mark Bulling.I was wondering > if you have discovered/developed a function to calculate the marginal effects > of the selection and outcome equations (sampleSelection). > > > > Thank you > > Marinella > [[alternative HTML version deleted]] > > __ > R-help@r-project.org 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] [External] Re: unable to access index for repository...
On Fri, 9 Oct 2020 at 12:28, Martin Maechler wrote: > > >>>>> Steven Yen > >>>>> on Fri, 9 Oct 2020 05:39:48 +0800 writes: > > > Oh Hi Arne, You may recall we visited with this before. I > > do not believe the problem is algorithm specific. The > > algorithms I use the most often are BFGS and BHHH (or > > maxBFGS and maxBHHH). For simple econometric models such > > as probit, Tobit, and evening sample selection models, old > > and new versions of R work equally well (I write my own > > programs and do not use ones from AER or > > sampleSekection). For more complicated models the newer R > > would converge with not-so-nice gradients while R-3.0.3 > > would still do nicely (good gradient). I use numerical > > graduent of course. I wonder whether numerical gradient > > routine were revised at the time of transition from > > R-3.0.3 to newer. > > As R-core member, particularly interested in numerical accuracy > etc, I'm also interested in learning what's going on here. > > I think we (R core) have never heard of anything numerically deteriorating > going from R 3.0.x to R 4.0.x, and now you are claiming that in > public, you should really post *reproducible* code giving > evidence to your claim. > > As was mentioned earlier, the difference may not be in R, but > rather in the versions of the (non-base R, but "extension") R > packages you use; and you were saying earlier you will check > that (using the old version of the 'maxLik' package with a newer > version of R and vice verso) and tell us about it. > > Thank you in advance on being careful and rational about such > findings. I totally agree with Martin: It would be good if Steven could run his code with different combinations of versions of R and maxLik so that we know whether the problem is caused by newer versions of R or by newer versions of maxLik (and hopefully also which version introduced the problem). Steven wrote that the problem appeared both in maxBFGS() and in maxBHHH(). This is somewhat surprising because maxBFGS() uses for the optimisation optim() which is implemented in R's base package "stats", while maxBHHH() uses for the optimisation maxNR() which is implemented in the maxLik package. Hence, it would be a very unlikely coincidence if the optimisation routines in optim() and maxNR() would become worse at the same time. I remember that we slightly changed the user interface and the default values of some arguments of some function in the maxLik package (e.g., maxBFGS, maxBHHH, maxNR, ...) a long time ago. I suggest that Steven checks whether he needs to update his code to reflect the changes in the user interface and/or in the default values of some arguments. Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] [External] Re: unable to access index for repository...
Hi Steven Which optimisation algorithms in maxLik work better under R-3.0.3 than under the current version of R? /Arne On Thu, 8 Oct 2020 at 21:05, Steven Yen wrote: > > Hmm. You raised an interesting point. Actually I am not having problems with > aod per se—-it is just a supporting package I need while using old R. The > essential package I need, maxLik, simply works better under R-3.0.3, for > reason I do not understand—specifically the numerical gradients of the > likelihood function are not evaluated as accurately in newer versions of R in > my experience, which is why I continue to use R-3.0.3. Because I use this > older version of R, naturally I need to install other supporting packages > such as aod and AER. > Certainly, I will install the zip file of the older version of maxLik to the > latest R and see what happens. Thank you. > > I will install the new maxLik in old R, and old maxLik in new R, and see what > happens. > > Sent from my iPhone > Beware: My autocorrect is crazy > > > On Oct 9, 2020, at 2:17 AM, Richard M. Heiberger wrote: > > > > I wonder if you are perhaps trying to solve the wrong problem. > > > > If you like what the older version of the aod package does, but not > > the current version, > > then I think the solution is to propose an option to the aod > > maintainer that would restore your > > preferred algorithm into the current version, and then use the current R. > > > > A less good, but possibly workable, option is to compile the old > > version of aod into the current R. > > > >> On Thu, Oct 8, 2020 at 1:45 PM Jeff Newmiller > >> wrote: > >> > >> All support on this list is voluntary, and support for old versions of R > >> is not even necessarily on-topic here which is why you keep getting nudged > >> to upgrade. Your "need" for support for an old version is definitely not > >> "our" problem, so I suggest you start looking for a consultant if this > >> issue is that important to you. Such is the nature of volunteer-developed > >> open source software... so support your local experts. > >> > >>> On October 8, 2020 10:22:54 AM PDT, Steven Yen wrote: > >>> Thanks for the help. I have a reason to continue with R-3.0.3. I used > >>> maxLik to estimate econometric models and some of them are better > >>> handled with R-3.0.3 (but not later)a sad reality I do not like. > >>> > >>> Here is what I did. I downloaded > >>> > >>> https://cran-archive.r-project.org/bin/windows/contrib/3.0/aod_1.3.zip > >>> > >>> and installed the zip file, which worked in both RStudio and R (without > >>> > >>> RStudio). > >>> > >>> In RStudio, I go Tools -> Install packages -> Install from -> (Choose > >>> zip) -> (Browse to the zip file) > >>> > >>> IN R, I go Packages -> Install packages from local file(s) -> (Browse > >>> to > >>> the zip file)... > > __ > R-help@r-project.org 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] Tobit Regression with unbalanced Panel Data
Dear Vanessa
Please provide a minimal *reproducible* example that illustrates your
problem, e.g. using a data set that is included in an R package.
Best regards,
Arne
On 4 January 2017 at 10:28, Vanessa Romero <vanrom...@gmail.com> wrote:
> Hello,
>
> I am doing Tobit Regression in R, because my dependent variable is censored
> at 0. I have unbalanced panel data, for 6 years, 107 companies. I use
> package CensReg.
>
> I have imported my database(T1).
>
> I use pdata.frame to specify the structure of my panel data. Like:
>
>
> *mydata<- pdata.frame (T1, index = c("firm", "year")) *
> Afterwards:
>
> *Tob <- censReg(formula=Imp ~ Bath + CEOTurnover + ChangeOCF + E + Sales +
> ROE + GTA + Size , data = mydata, method="BHHH") *
> (as explained here:
> https://cran.r-project.org/web/packages/censReg/vignettes/censReg.pdf)
>
> I got here error message:
>
>
> *Warnmeldung: In log(rEff$ercomp$sigma$id) : NaNs wurden erzeugt*
>
> Another error message when *summary(Tob)*
>
>
>
>
>
> *Call: censReg(formula = Imp ~ Bath + CEOTurnover + ChangeOCF + E + Sales +
> ROE + GTA + Size, data = mydata, method = "BHHH") Observations: Total
> Left-censored Uncensored Right-censored 606 469 137 0 Coefficients: Fehler
> in printCoefmat(coef(x, logSigma = logSigma), digits = digits) : 'x' must
> be coefficient matrix/data frame*
>
> I am new to statistics and to R, what could be the problem or would you
> suggest using other package.
>
> Thank you,
> Vanessa
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org 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.
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org 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.
Re: [R] Hickman models with two binary dependent variables in R
Dear Faradj On 27 August 2016 at 16:11, Faradj Koliev <farad...@gmail.com> wrote: > It actually worked with heckit() command as well, do I need to use > selection()? I suggest that you use the selection() command/function only. heckit() is just a wrapper to selection() with the only difference that heckit() estimates the model with the 2-step method by default, while selection() estimates the model by the maximum likelihood method by default (unless specified otherwise with argument "method"). I see no good reason for estimating the sample-selection model by the 2-step method instead of the maximum likelihood method. > Also, I would be really grateful if you can suggest a package that would > allow for estimation of heckman models with two ordered variables (0-1-2). > Can sampleSelection handle this? Do you mean that the dependent variable of the selection equation is an ordered variable with three levels (0-1-2), that the dependent variable of the outcome equation is an ordered variable with three levels (0-1-2), or that both of the dependent variables are ordered variable with three levels (0-1-2) each? In any case: no, this is not implemented in the sampleSelection package. Sorry! Anyway, I guess that it is not too complicated to derive the likelihood function and implement the estimation yourself, e.g., using the maxLik package. If you do this, I would be happy to help you to implement this feature in the sampleSelection package. Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] Hickman models with two binary dependent variables in R
ing any reasoning. > - if the e-mail contains an offer, the recipient is entitled to immediately > accept such offer; The sender of this e-mail (offer) excludes any acceptance > of the offer on the part of the recipient containing any amendment or > variation. > - the sender insists on that the respective contract is concluded only upon > an express mutual agreement on all its aspects. > - the sender of this e-mail informs that he/she is not authorized to enter > into any contracts on behalf of the company except for cases in which he/she > is expressly authorized to do so in writing, and such authorization or power > of attorney is submitted to the recipient or the person represented by the > recipient, or the existence of such authorization is known to the recipient > of the person represented by the recipient. > __ > R-help@r-project.org 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] HELP - Hausman Test + systemfit
> mention in *"inst"*. I put additional variable in Model2 which does not >>> exist in Model1, but I ma not sure that is correct! Or maybe I need to >>> put >>> all variables, which were used in 2 models. if You can explain - thank >>> you! >>> system <- list(Model1, Model2) >>> >>> # perform the estimations >>> fit2sls <- systemfit(system, "2SLS", inst = inst, data = ACU) >>> >>> but R responded: >>> >>> Error in systemfit(system, "2SLS", inst = inst, data = ACU) : >>> the list of argument 'formula' must contain only objects of class >>> 'formula' >>> >>> >>> >>> Please, help me to understand What I do wrong! >>> Best, >>> Kateryna >>> >>> [[alternative HTML version deleted]] >>> >>> __ >>> R-help@r-project.org 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]] >> >> __ >> R-help@r-project.org 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. > > David Winsemius > Alameda, CA, USA > > __ > R-help@r-project.org 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] Estimating Endogenous Selection Model
Dear Johannes
You can use function selection() of the sampleSelection package and
estimate an endogenous switching regression ("tobit-5") model, where
you have the selection equation to model investmentment vs.
non-investment and you have two outcome equations, one for firms that
invest and the other for firms that do not invest.
Best regards,
Arne
On 18 September 2015 at 07:43, Johannes Muck <m...@econ.uni-frankfurt.de> wrote:
> Dear all,
>
> I want to estimate a model in which individuals self-select into two
> different actions (e.g. invest or not invest). Moreover, the factors that
> influence the selection decision also affect the ultimate outcome variable
> (e.g. return on investment). That is, I want to estimate a model with
> endogenous selection.
>
> My question is: Which R package supports this kind of estimation?
> To me it seems like the package sampleSelection
> (https://cran.r-project.org/web/packages/sampleSelection/sampleSelection.pdf
> ) can only be used in case of a SAMPLE selection, i.e. in the example above:
> only individuals that invest are observed. Yet in my case, both actions are
> observed but the decision which action to choose is endogenous. Or does this
> package also support self-selection models?
>
> Thank you very much in advance.
>
> Best,
>
> Johannes
>
> __
> R-help@r-project.org 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.
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org 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.
Re: [R] simultaneous equation model with endogenous interaction terms
Dear Janka On 10 August 2015 at 11:25, Janka Vanschoenwinkel janka.vanschoenwin...@uhasselt.be wrote: Dear list members, I am building a model such as: Y1 = Y2*X1 + X2 Y2 = Y1*X1 + X2 Do you mean the model: Y1 = b10 + b11 * (Y2*X1) + b12 * X2 + e1 Y2 = b20 + b21 * (Y1*X1) + b22 * X2 + e2 where Y1 and Y2 are two (endogenous) dependent variables, X1 is a potentially endogenous explanatory variable, X2 is an exogenous explanatory variable, e1 and e2 are two potentially contemporaneously correlated error terms, and b10, b11, b12, b20, b21, and b22 are parameters to be estimated? X2 is the exogenous variable Z1 is the instrument of Y1 Z2 is the instrument of Y2 This is a simultaneous equation model. I know how to build a simultaneous equation model without interaction terms: library(systemfit) eq1 - Y1~Y2+X2+Z2 eq2 - Y2~Y1+X2+Z1 inst - ~X2+Z1+Z2 system - list(eq1=eq1, eq2=eq2) reg2SLS -systemfit(system, 2SLS, inst=inst, data=mydata) summary(reg2SLS) I also know how to do a normal 2SLS with interaction terms: library(systemfit) ivreg(Y1~Y2*X1 | Z2*X1, data= Alldata) However, I don't know how to deal with the interaction terms in the simultaneous equation model. I am experimenting both with R and STATA to see which formulation gives the same result in both softwares, but until know without success. Could somebody help me with this? To estimate the above model specification, the following should work: eq1 - Y1 ~ I(Y2*X1) + X2 eq2 - Y2 ~ I(Y1*X1) + X2 inst - ~ X2 + Z1 + Z2 system - list( eq1 = eq1, eq2 = eq2 ) reg2SLS - systemfit( system, 2SLS, inst = inst, data = mydata ) Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] Warning message with maxLik()
the comparederivitive() function, and the analytic and numeric gradients were so close. Any help please? Maram Salem [[alternative HTML version deleted]] __ R-help@r-project.org 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] NaN produced from log() with positive input
Dear Maram On 8 July 2015 at 17:52, Maram Salem marammagdysa...@gmx.com wrote: Dear Arne, On a second thought, as per your mail the warning messages occur each time, when maxLik() tries to calculate the logLik value for theta[1] = 0, theta[1] + theta[2] = 0, theta[3] = 0 or something similar. The component of the theta vector are all indeed strictly positive, and the initial values I used are c(40,50,2). and this means that theta[1] 0, theta[1] + theta[2] 0, and theta[3] 0 . These initial values are the parameter values that were used for generating the data (the C and T vectors). That's why I don't know why the warnings occur in the first place and why the estimates are far away from the initial values.a Any suggestions? - don't send your messages twice; - do what I suggested in my previous e-mail; - increase the number of observations; - check your log-likelihood function; - use an optimisation algorithm that does not use derivatives; - use numeric derivatives in the optimisation; - check your function for returning the derivatives of the log-likelihood function, e.g. with compareDerivatives(); - check the data generating process Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] NaN produced from log() with positive input
summary (mle), I got Maximum Likelihood estimation Newton-Raphson maximisation, 7 iterations Return code 1: gradient close to zero Log-Likelihood: -55.89012 3 free parameters Estimates: Estimate Std. error t value Pr( t) [1,] 11.132Inf 0 1 [2,] 47.618Inf 0 1 [3,]1.293Inf 0 1 Where the estimates are far away from the starting values and they have infinite standard errors. I think there is a problem with my gradlik or hesslik functions, but I can't figure it out. Any help? Thank you in advance. Maram [[alternative HTML version deleted]] __ R-help@r-project.org 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] PEA and APE Tobit
Dear Annelies On 26 February 2015 at 09:12, hnlki annelies.hoebe...@ugent.be wrote: I estimated a tobit model tobit.fit-tobit(y~x,left=0, right=Inf) (library AER) or tobit2.fit-censReg(y~x, left=0, right=Inf) (librarycensReg) I' have estimated the partial effect at the average as: pea-(pnorm((colMeans(x)%*%tobit.fit$coef[-1])/tobit.fit$scale))%*%tobit.fitt$coef[-1] and the average partial effect as: ape- (length(x[,1]))^(-1)*sum(pnorm((x%*%tobit.fit$coef[-1])/tobit.fit$scale))*tobit.fit$coef[-1] I guess I did something wrong as margEff( tobit2.fit) (library(censReg) gives a different result than my partial effect at the average. Any ideas about what I did wrong? I did not find the underlying code of margEff. [...] PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Please follow the posting guide and provide a self-contained reproducible example. Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] MLE
Dear Pari On 21 December 2014 at 06:59, pari hesabi statistic...@hotmail.com wrote: I am trying to get the Maximum Likelihood Estimation of a parameter in a probability mass function. The problem is my pmf which includes a summation and one integral. it is not similar to other known pdfs or pmfs such as normal, exponential, poisson, . Does anybody know whether I can use the current packages(like Maxlik) in R for getting the MLE of the parameter? maxLik (not Maxlik) will probably work. Can anybody explains me with an example? R library( maxLik ) R ?maxLik http://dx.doi.org/10.1007/s00180-010-0217-1 https://absalon.itslearning.com/data/ku/103018/publications/maxLik.pdf I would appreciate any help. http://www.R-project.org/posting-guide.html http://maxlik.org/ https://r-forge.r-project.org/projects/maxlik/ Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org 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.
Re: [R] maximum likelihood estimation
On 10 October 2014 08:04, pari hesabi statistic...@hotmail.com wrote:
Hello,As an example for Exponential distribution the MLE is got by this
structure:t - rexp(100, 2)loglik - function(theta){ log(theta) - theta*t}a
- maxLik(loglik, start=1)print(a)Exponential distribution has a simple
loglikelihood function. But if a new pdf has a more complicated form
like:pr(N(2)=n)= integral( ((2-x)^n)*(exp(ax-2))) - integral (((5-ax)^n)),
both integrals are defined over the interval(0,2) with respect to x. I
also need to use the loglike of the form : [F log(pr(n))]=lnL where F is the
vector of observations and (n) is the vector of input for the defined pmf.
Do you think how I can insert (define) the pr(n) in the loop below? loglik
- function(a) sum(F*(log(pr(n)))? n - c(0,1,2,3,4,5,6,7,8)
F-c(0,0,1,3,5,7,8,11,10) loglik - function(a) sum(F*(log(pr(n)))??
re - maxLik (loglik, start=.5) summary(re)I would be grateful if you
let me know your idea.Best Regards,pari
[...]
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Your example is not reproducible:
pr(N(2)=n)= integral( ((2-x)^n)*(exp(ax-2))) - integral (((5-ax)^n))
Error: unexpected '=' in pr(N(2)=
Please format your email in a way that one can easily copy the R code
to the R console.
Best regards,
Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] maximum likelihood estimation
Dear Pari
On 7 October 2014 10:55, pari hesabi statistic...@hotmail.com wrote:
HelloI am trying to estimate the parameter of a function by the Maximum
Likelihood Estimation method.If the function is the difference between two
integrals: C-function(n){integrand3-function(x) {((2-x)^n)*(exp(ax-2))}cc-
integrate (integrand3,0,2)print(cc)}
D-function(n){integrand4-function(x) {((5-ax)^n)}cc- integrate
(integrand4,0,2)print(cc)}
f(n) = C(n) - D(n)
I need to estimate parameter (a). loglikelihood function is in the form of
sum[F(k) log(f(k))]=lnL. I am wondering how to introduce my logliklihood
function to the loop below. Can anybody help me for correcting the following
loop? if there are some other packages better than this , please let me
know.Thank you,Diba
n - c(0,1,2,3,4,5,6,7,8)
F-c(0,0,1,3,5,7,8,11,10)
loglik - function(a) sum(F*log(C(n)-D(n)))
re - maxLik (loglik, start=.5)
summary(re)
Unfortunately, I cannot reproduce your example, because I get an error
message about an unexpected symbol.
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Best regards,
Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] comparing two half-normal production stochastic frontier functions
Dear Rainer On 3 October 2014 14:51, Rainer M Krug rai...@krugs.de wrote: I am using the function frontier::sfa (from the package frontier) to estimate several half-normal production stochastic frontier functions. Now I want to compare the coefficients of the linear frontier function and see if they are different. According to my stackexchange (CrossValidated) question [1] I can compare these as I can compare a normal linear regression. In R, I would uswe the function anova to do this model comparison - correct? Now this function does not accept objects of the type 'frontier' - so how can I do this comparison in R? To re-iterate, I want to know if the coefficients of the frontier line (slope and intercept) are significantly different. Below please find a reproducible example based on data provided in the package, of what I did, and below the transcript. --8---cut here---start-8--- library(frontier) data(front41Data) dat1 - front41Data[1:30,] dat2 - front41Data[30:60,] x1 - sfa(log(output) ~ log(capital), data=dat1) x2 - sfa(log(output) ~ log(capital), data=dat2) x1 x2 anova(x1, x2 --8---cut here---end---8--- library( frontier ) data( front41Data ) # estimate pooled model mp - sfa( log(output) ~ log(capital), data = front41Data ) # create a dummy variable front41Data$dum - rep( c( 1, 0 ), 30 ) # estimate model with different intercepts and different slopes # but the same sigmsSq and the same gamma md - sfa( log(output) ~ log(capital)*dum, data = front41Data ) # likelihood ratio test lrtest( mp, md ) If you have further questions regarding the frontier package, you may also use the help forum at frontier's R-Forge site: https://r-forge.r-project.org/projects/frontier/ ...and please do not forget to cite the frontier package in your publications (see output of the R command 'citation(frontier)'). Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
[R] [R-pkgs] sampleSelection 1.0 released
Dear All The recently released version of the sampleSelection package [1] has the version number 1.0, because we consider the package now to be mature and feature-complete. Compared to version 0.6 of the sampleSelection package (which is described in an article in the Journal of Statistical Software [2]), the new version has, for instance, following new features: * estimation of binary-outcome selection models * probit models, tobit-2 models, and binary-outcome selection models can now be estimated with observation-specific weights * predict() methods for calculating unconditional and conditional expectations * improved/extended residuals() methods * logLik() methods * nobs() and nObs() methods * default df.residual() method works * removed the deprecated function tobit2() Any kind of feedback is--as always--very welcome, preferable via R-Forge [3]! Best regards, Ott and Arne [1] http://cran.r-project.org/package=sampleSelection [2] http://www.jstatsoft.org/v27/i07/ [3] http://r-forge.r-project.org/projects/sampleselection/ -- Arne Henningsen http://www.arne-henningsen.name ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@r-project.org mailing list 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.
Re: [R] Tobit model with panel data
Hi Phil On 30 April 2014 23:06, phil quij...@gmx.net wrote: I just started to work with R a couple of weeks ago. Right now I would like to regress an independent variable on a couple of explanatory variables. The dependent variable is left censored in the sense that all negative values and zero are set equal to one. That is most likely not a suitable data transformation. This is done because I want to take the logarithm which, as you all know, is only defined for values bigger than zero. However, what I don't know now is how to compute this on R. Does anybody know how to proceed? I already checked http://cran.r-project.org/web/packages/censReg/vignettes/censReg.pdf http://cran.r-project.org/web/packages/censReg/vignettes/censReg.pdf but I am not able to apply the commands shown on page 8 correctly on my data set. Hence, I would really appreciate if somebody could give me a step-by-step instruction for my own dataset. 1st step: load the data set into R (see http://cran.r-project.org/doc/manuals/r-release/R-data.html) 2nd step: use pdata.frame() or plm.data() to specify the panel structure (see http://www.jstatsoft.org/v27/i02/) 3rd step: use censReg() to estimate the model (see http://cran.r-project.org/web/packages/censReg/vignettes/censReg.pdf) Please do read the posting guide http://www.R-project.org/posting-guide.html Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Problem with new(er) R version's matrix package
On 25 April 2014 20:15, David Winsemius dwinsem...@comcast.net wrote: On Apr 25, 2014, at 9:17 AM, Werner W. wrote: Dear Rs, I am re-executing some older code. It does work in the ancient R 2.12.0 which I still have on my PC but with the new version R 3.1.0 it does not work any more (but some other new stuff, which won't work with 2.12). The problem arises in context with the systemfit package using the matrix package. In R 3.1.0 the following error is thrown: Error in as.matrix(solve(W, tol = solvetol)[1:ncol(xMat), 1:ncol(xMat)]) : error in evaluating the argument 'x' in selecting a method for function 'as.matrix': Error in .solve.sparse.dgC(as(a, dgCMatrix), b = b, tol = tol) : LU computationally singular: ratio of extreme entries in |diag(U)| = 7.012e-39 However, I have no clue what I can do about this. Was there some change in the defaults of the matrix package? I couldn't find anything apparent in the changelog. As the same code works in R 2.12.0, I suppose that the problem is not my data. You have not told us what version of the Matrix package you were using. As such I would suggest that you review the Changelog which is a link for the CRAN page for pkg:Matrix and go back 4 years or so since R major versions change about once a year. http://cran.r-project.org/web/packages/Matrix/ChangeLog In addition, please provide a minimal, self-contained, reproducible example. Best, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Testing correlation of equation in a SUR model fitted by systemfit
Dear Paul On 15 April 2014 19:23, Paul Smith phh...@gmail.com wrote: How to test whether the correlation in the matrix of correlation of a two-equations SUR model fitted by package systemfit are significant? You can use a likelihood-ratio test to compare the SUR model with the corresponding OLS model. The only difference is that the OLS model assumes that all off-diagonal elements of the covariance matrix of the residuals of the different equations are zero. An example: library( systemfit ) # load data set data( Kmenta ) # specify system of equations eqDemand - consump ~ price + income eqSupply - consump ~ price + farmPrice + trend system - list( demand = eqDemand, supply = eqSupply ) # estimate OLS model fitols - systemfit( system, data=Kmenta ) # estimate SUR model fitsur - systemfit( system, SUR, data = Kmenta ) # LR test: OLS vs. SUR lrtest( fitols, fitsur ) # estimate iterated SUR model fititsur - systemfit( system, SUR, data = Kmenta, maxit = 100 ) # LR test: OLS vs. SUR lrtest( fitols, fititsur ) If you have further questions regarding the systemfit package, you can also use the help forum at systemfit's R-Forge site: https://r-forge.r-project.org/projects/systemfit/ ... and please do not forget to cite systemfit in your publications (see output of the R command 'citation(systemfit)'). Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Discrete-continuous equation system in R
Dear Reinhard
On 2 January 2014 14:12, Reinhard Hössinger
reinhard.hoessin...@boku.ac.at wrote:
I want to estimate an equation system with 3 nonlinear continuous
equations and one discrete choice (using multinomial logit for the
latter). For the nonlinear continuous equations, function nlsystemfit
{systemfit} seems to be appropriate. But what's about the discrete
choice? Its error component has a logistic distribution. Can it still be
incorporated in the equation system? Or can/must the error component be
treated in some way to be incorporated? I cannot find a reference to
discrete choices in the systemfit package description nor somewhere
else.
I guess that the econometric estimation of this type of model has not
been implemented in any ready-to-use software package. Hence, you
probably need to derive the log-likelihood function of this model,
implement this function (and preferably also the derivatives/gradients
with respect to the parameters), and maximize this log-likelihood
function, e.g. with the R package maxLik (http://maxlik.org/).
Best regards,
Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] MaxLik estimation issues
Dear Filipe
On 14 October 2013 17:42, Filipe Ribeiro flipjribe...@hotmail.com wrote:
Dear Arne,
First of all, thank you very much for your reply. Secondly, please accept my
apologies for only give you an answer now, but I was moving from one country
to another and only today I was able to get back to work.
I did not write the model specification because it was just a trial, but my
main idea was to apply the Nelder-Mead model.
Since the beginning my guess was that something is wrong with the
log-likelihood function, however, I don't find any problem with the same
exact function applying the optim function:
loglike.GGompi - function(theta,age,deaths,exposures) {
alpha - exp(theta[1])
beta - exp(theta[2])
gamma - exp(theta[3])
first - alpha*exp(beta*age)
second - 1+(((alpha*gamma)/beta)*(exp(beta*age)-1))
mu - first/second
llk - -sum((deaths * log(mu)) + (- mu*exposures))
return(llk)
}
fit1 - optim(par=c(-4.1402817, -0.6375773, -1.6945914),
loglike.GGompi,
age=0:6,
deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
231951.64, 69502.65, 15798.72))
exp(fit1$par)
exp(fit1$par)
[1] 0.01585683 0.53471945 0.25368426
Due to this results I can't understand why...
Please note that optim() is (by default) *minimizing* the provided
function, while maxLik() is *maximizing* the provided log-likelihood
function.
Cheers,
Arne
A 26/09/2013, às 05:28, Arne Henningsen escreveu:
Dear Filipe
On 25 September 2013 14:23, Filipe Ribeiro flipjribe...@hotmail.com wrote:
Hello everybody!
I'm having some trouble to compute maximum likelihood
estimations using maxLik package and I hope that you
could give me a hint.
The main problem is that I'm not able to get a result not
even close to the ones given by glm() directly, and the
second one is: Error in maxNRCompute(fn = logLikAttr,
fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in
gradient.
The codes:
loglike.GGompiMaxLik - function(theta,age,deaths,exposures) {
alpha - exp(theta[1])
beta - exp(theta[2])
gamma - exp(theta[3])
first - alpha*exp(beta*age)
second - 1+(((alpha*gamma)/beta)*(exp(beta*age)-1))
mu - first/second
llk - -sum((deaths * log(mu)) + (- mu*exposures))
return(llk)
}
fit1 - maxLik(loglike.GGompiMaxLik,
age=0:6,
deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
231951.64, 69502.65, 15798.72),
start=c(-4.1402817, -0.6375773, -1.6945914))
Do you know how I can solve this problem?
You did not write which model specification you want to estimate but I
am pretty sure that something in your log-likelihood function is
incorrect. The log-likelihood value at the starting values of the
parameters is so large that R even cannot calculate the likelihood
value:
a - loglike.GGompiMaxLik(c(-4.1402817, -0.6375773, -1.6945914), age=0:6,
+ deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
+ exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
+ 231951.64, 69502.65, 15798.72))
a
[1] 580365.2
exp(a)
[1] Inf
In the second iteration, the first parameter gets so small (large in
absolute terms, -5e+10) that the log-likelihood value become extremely
(numerically infinitely) large and the gradients cannot be computed
(by the finite-difference method):
fit1 - maxLik(loglike.GGompiMaxLik,
+ age=0:6,
+ deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
+ exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
+ 231951.64, 69502.65, 15798.72),
+ start=c(-4.1402817, -0.6375773, -1.6945914))
Iteration 2
Parameter:
[1] -5.174233e+10 -3.839076e+02 5.988668e+00
Gradient:
[,1] [,2] [,3]
[1,] NaN NaN NaN
Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad,
hessOrig = hess, :
NA in gradient
b - loglike.GGompiMaxLik(c(-5.174233e+10, -3.839076e+02, 5.988668e+00),
age=0:6,
+ deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
+ exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
+ 231951.64, 69502.65, 15798.72))
b
[1] Inf
Please note that you can also find hints and ask questions about the
maxLik package in the forums at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
...and please do not forget to cite the maxLik package in your publications:
http://cran.r-project.org/web/packages/maxLik/citation.html
Best wishes,
Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] MaxLik estimation issues
Dear Filipe
On 25 September 2013 14:23, Filipe Ribeiro flipjribe...@hotmail.com wrote:
Hello everybody!
I'm having some trouble to compute maximum likelihood
estimations using maxLik package and I hope that you
could give me a hint.
The main problem is that I'm not able to get a result not
even close to the ones given by glm() directly, and the
second one is: Error in maxNRCompute(fn = logLikAttr,
fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in
gradient.
The codes:
loglike.GGompiMaxLik - function(theta,age,deaths,exposures) {
alpha - exp(theta[1])
beta - exp(theta[2])
gamma - exp(theta[3])
first - alpha*exp(beta*age)
second - 1+(((alpha*gamma)/beta)*(exp(beta*age)-1))
mu - first/second
llk - -sum((deaths * log(mu)) + (- mu*exposures))
return(llk)
}
fit1 - maxLik(loglike.GGompiMaxLik,
age=0:6,
deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
231951.64, 69502.65, 15798.72),
start=c(-4.1402817, -0.6375773, -1.6945914))
Do you know how I can solve this problem?
You did not write which model specification you want to estimate but I
am pretty sure that something in your log-likelihood function is
incorrect. The log-likelihood value at the starting values of the
parameters is so large that R even cannot calculate the likelihood
value:
a - loglike.GGompiMaxLik(c(-4.1402817, -0.6375773, -1.6945914), age=0:6,
+ deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
+ exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
+ 231951.64, 69502.65, 15798.72))
a
[1] 580365.2
exp(a)
[1] Inf
In the second iteration, the first parameter gets so small (large in
absolute terms, -5e+10) that the log-likelihood value become extremely
(numerically infinitely) large and the gradients cannot be computed
(by the finite-difference method):
fit1 - maxLik(loglike.GGompiMaxLik,
+ age=0:6,
+ deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
+ exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
+ 231951.64, 69502.65, 15798.72),
+ start=c(-4.1402817, -0.6375773, -1.6945914))
Iteration 2
Parameter:
[1] -5.174233e+10 -3.839076e+02 5.988668e+00
Gradient:
[,1] [,2] [,3]
[1,] NaN NaN NaN
Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad,
hessOrig = hess, :
NA in gradient
b - loglike.GGompiMaxLik(c(-5.174233e+10, -3.839076e+02, 5.988668e+00),
age=0:6,
+ deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198),
+ exposures=c(935023.67, 819392.00, 724568.17, 470947.00,
+ 231951.64, 69502.65, 15798.72))
b
[1] Inf
Please note that you can also find hints and ask questions about the
maxLik package in the forums at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
...and please do not forget to cite the maxLik package in your publications:
http://cran.r-project.org/web/packages/maxLik/citation.html
Best wishes,
Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] Weighted SUR/NSUR
Dear Ariel Thank you for your detailed explanations and the example. Indeed, it should be rather straightforward to implement observation-specific weights in systemfit (i.e. in the estimation of systems of linear or non-linear equations). As you indicated that many people are looking for this feature, I wonder if I should start a crowd-funding campaign... Best regards, Arne On 16 August 2013 18:23, Ariel ariel.muld...@oregonstate.edu wrote: Arne Henningsen-3 wrote Is it possible to run SUR with weights using systemfit? I mean weighted seemingly unrelated regression (weighted SUR) Currently, systemfit cannot estimate (SUR) models with observation-specific weights :-( or weighted nonlinear unrelated regression (weighted NSUR). We are still not yet finished with implementing nonlinear models in systemfit (see http://www.systemfit.org/) :-( I recently had a student come to me with a very similar (okay, identical) problem as the OP. I had to learn PROC MODEL, anyway, so I thought I’d poke around in R while I was at it. I have nothing to add about any problems with or the lack of maturity of the estimation procedure for nlsystemfit(), but I do have some ideas about observation-level weights. It took me awhile to make the leap from the fairly straightforward linear weighted least squares (for example, see Weisberg's Applied Linear Regression textbook equation 5.8) to understanding how weighting worked in nonlinear least squares. The R help forum certainly came in handy: https://stat.ethz.ch/pipermail/r-help/2004-November/060424.html. I can add weights into a nonlinear regression by simply multiplying both the response and the nonlinear function by the square root of the desired weights. Here’s a toy example, where I compare a model fit using the “weights” argument in nls() with a model where I put the weights in “by hand” : DNase1 = subset(DNase, Run == 1) fit2 = nls(density ~ Asym/(1 + exp((xmid - log(conc))/scal)), data = DNase1, start = list(Asym = 3, xmid = 0, scal = 1), weights = rep(1:8, each = 2)) summary(fit2) # Take the square root of the weights for fitting “by hand” sw = sqrt(rep(1:8, each = 2) ) fit3 = nls(sw*density ~ sw*(Asym/(1 + exp((xmid - log(conc))/scal))), DNase1, start = list(Asym = 3, xmid = 0, scal = 1) ) summary(fit3) # The predicted values for fit3 need to be divided by the weights # but the residuals are weighted residuals predict(fit2) predict(fit3)/sw It seems like this weighted approach could be easily extended to the model formulas for a system of nonlinear equations (it would be similar for linear equations) to be fit with systemfit. Parresol (2001) in his paper titled Additivity of nonlinear biomass equations has run weighted NSUR using PROC MODEL (SAS institute Inc.1993). I was wondering if r can do that. It turned out I had to use this weighting approach in PROC MODEL, as well, when each equation in the system had a different set of weights. The estimates I get when fitting the Parresol example mentioned by the OP using nlsystemfit and PROC MODEL are within spitting distance of each other, so I feel like I am at least making the same mistakes in both software packages. I'm wondering if my logic is sound or if I'm missing some complication that occurs when working with systems of equations. I’ve seen several folks looking to fit weighted systems of equations in R with systemfit, and this approach might get them what they need. Ariel -- View this message in context: http://r.789695.n4.nabble.com/Weighted-SUR-NSUR-tp4670602p4673973.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] 2sls and systemfit
Dear Cecilla! On 22 July 2013 14:08, Cecilia Carmo cecilia.ca...@ua.pt wrote: I have the following model: Cost of debt = intercept + information quality + control variable1 + control variable2 + … + error term I want to perform 2sls because I think I could have two situations: First: maybe information quality is correlated with the error term (because of omitted variables) Second: maybe information quality depends on the cost of debt like this: information quality = intercept + Cost of debt + control variable1 + control variable2 + … + error term I have some variables (instruments) for information quality. Now I need to know how to use systemfit with this information in each of those situations. The use of systemfit depends on the econometric model that you want to estimate. In the first case, the following specification might be suitable: result1 - systemfit( Cost_of_debt ~ intercept + information_quality + control_variable1 + control_variable2, inst = ~ control_variable1 + control_variable2 + instrument_variable1 + instrument_variable2, method = 2SLS ) Please note that you have to assume that the instrumental variables (instrument_variable1 and instrument_variable2) do not influence the cost of debt. In the second case, you can only estimate the model if you have some exclusion restrictions, i.e. at least one exogenous variable that does not influence the cost of debt (e.g. instrument_variable1 and instrument_variable2) and at least one exogenous variable that does not influence the information quality (e.g. control_variable2). In this case, the following specification might be suitable: sys - list( Cost_of_debt ~ intercept + information_quality + control_variable1 + control_variable2, information_quality ~ intercept + Cost_of_debt + control_variable1 + instrument_variable1 + instrument_variable2 ) result2 - systemfit( sys, inst = ~ control_variable1 + control_variable2 + instrument_variable1 + instrument_variable2, method = 2SLS ) You could also use the 3SLS method in this case. Please do not forget to cite the systemfit package in your publications. Thanks! Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
[R] [R-pkgs] Version 1.0 of the R package frontier released
Dear all I am happy to announce that version 1.0 of the frontier package is available on CRAN. The R package frontier provides tools for analysing efficiency and productivity using the stochastic frontier approach. This R package is based on Tim Coelli's DOS software FRONTIER 4.1 and has been available on CRAN for almost 5 years now. After many improvements of the source code, a few additional features, and a few changes of the user interface in the beginning, the code and the user interface have been rather stable in the previous 2.5 years. The frontier package has been used by many R users and many applications have proven its reliability. Therefore, I have called the latest version of this package 1.0. This version is almost identical to the previous versions but it includes citation information and the argument farrell of efficiencies.frontier() has been renamed as minusU and the argument farrell of summary.frontier() has been renamed as effMinusU, because the term farrell was misleading in some cases. The argument farrell can still be used for maintaining backward-compatibility. http://cran.r-project.org/package=frontier https://r-forge.r-project.org/projects/frontier/ Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@r-project.org mailing list 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.
Re: [R] 2SLS / TSLS / SEM non-linear
Dear HC On 30 June 2013 18:53, hck hans-christian.krumh...@uni-ulm.de wrote: Generally I would have the following equations X_i = IV3_i + IV4_i * Y_i applying for every company (i). In a first step, I am interested in estimating the relationship between X and Y: Y_i = a + b * X_i + u to ultimatly estimate X_i by substituting the Y_i and solving for X_i to be able to estimate the X_i by just IV3_i, IV4_i, and the a and b. Now, let's construct values from a sample of listed companies. In the capital market, I can observe IV3_i, IV4_i, and X_i. With these I calculate Y_i: Y_i = IV1_i + IV2_i * X_i (note: IV3 and IV4 are just a transformation of IV1 and IV2). Of course, I could rewrite this equation as Y_i = c + d * IV1_i + e * IV2_i * X_i + v. For a couple of observations, I have now combinations of X_i and Y_i to get the a and b coefficient by estimating Y_i = a + b * X_i + u. It seems to me that this estimation is very simple: myModel - lm( Y ~ X ) but perhaps I did not completely understand your model specification. Best, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Testing for weak exogeneity in a SUR ECM
Dear Kathrinchen It seems to me that your question is about statistics rather than about R and systemfit. If you find out how the statistical test should be conducted theoretically, I can probably advise you how to implement the test in R (systemfit). Best wishes, Arne On 11 July 2013 13:21, Kathrinchen katha.m...@web.de wrote: Dear all, I have set up a Labour Demand Error Correction Model for some German federal states. As I expect the labour markets to be correlated I used a Seemingly Unrelated Regression using systemfit in R. My Model is: d(emp)_it = c + alpha*ln(emp)_i,t-1 + beta_1*ln(gdp)_i,t-1 + + beta_2*ln(wage)_i,t-1 + + beta_1*ln(i)_i,t-1 + gamma_1*d(gdp)_it + gamma_2*d(wage)_it with emp_it being the employment in state i at time t, i stands for the real interest rate, ln() is the logarithmed data, while d() stands for the difference operator. I would like to test now for weak exogeneity and I am not quite sure what kind of regression to run. If I run: d(gdp)_it = c + alpha*ln(emp)_i,t-1 + beta_1*ln(gdp)_i,t-1 + + beta_2*ln(wage)_i,t-1 + + beta_1*ln(i)_i,t-1 + gamma_1*d(emp)_it + gamma_2*d(wage)_it with Systemfit, alpha is statistically significant, so I have to reject the hypothesis of weak exogeneity...Literature is in my opinion not so clear on what to test! I use data from an application, they conclude that endogeneity is not a problem: they regress the possible endogenous variables on the presumed equilibrium relation, a constant and one autoregressive lag - here I am not sure, what they mean. I would very much appreciate your help! Thanks a lot! -- View this message in context: http://r.789695.n4.nabble.com/Testing-for-weak-exogeneity-in-a-SUR-ECM-tp4671321.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] 2SLS / TSLS / SEM non-linear
Dear HC On 30 Jun 2013 13:02, hck hans-christian.krumh...@uni-ulm.de wrote: I started having a look into the systemfit. Basically, my two equations are: 1. Y = IV1 + IV2 x X 2. Y = a + b x X + u where Y and X are the two endogenous variables, IV1 and IV2 are the instruments, and u is the error term. I do not understand your model specification. In your first equation IV1 and IV2 look like parameters. Is your model perhaps: Y = a + b x X + u X = c + d x IV1 + e x IV2 + v with a, b, c, d, and e being parameters and u and v being disturbance terms? Best wishes, Arne [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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.
Re: [R] Weighted SUR/NSUR
Dear Tarquinio On 29 Jun 2013 21:01, Tarquinio magalhaes tarq...@yahoo.com.br wrote: Is it possible to run SUR with weights using systemfit? I mean weighted seemingly unrelated regression (weighted SUR) Currently, systemfit cannot estimate (SUR) models with observation-specific weights :-( or weighted nonlinear unrelated regression (weighted NSUR). We are still not yet finished with implementing nonlinear models in systemfit (see http://www.systemfit.org/) :-( Parresol (2001) in his paper titled Additivity of nonlinear biomass equations has run weighted NSUR using PROC MODEL (SAS institute Inc.1993). I was wondering if r can do that. Unfortunately not (see above). As several people have asked for observation-specific weights and NSUR, I should perhaps start a crowd-funding initiative for implementing these features. BTW: You can ask questions that are specific to systemfit via the help forum at systemfit's R-Forge site: https://r-forge.r-project.org/projects/systemfit/ Best wishes, Arne [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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.
Re: [R] maxLik: different estimations from different packages
Dear Alfonso
On 10 May 2013 12:51, alfonso.carf...@uniparthenope.it wrote:
we are computing maximum likelihood estimations using maxLik package and we
realized that the results of the estimation depend on the version of the
package installed
for example if we try to estimate this function with an old version of
maxLik under R 2.13.1 (32 bit version installed 2 years ago):
L-function (param) {b0t-param[1]
p1t-param[2]
p2t-param[3]
p3t-param[4]
p4t-param[5]
for(i in 17:T) {n[i,]- b0t + p1t*a[i-1] + p2t*sum(a[(i-4):(i-1)]) +
p3t*(sum(a[(i-8):(i-1)])) + p4t*(sum(a[(i-16):(i-1)]))
m[i,]-exp(n[i])/(1+exp(n[i]))
ll[i-16,]-a[i]*log(m[i])+(1-a[i])*log(1-m[i]) }
sum(ll)}
b2-maxLik(L, start=c(-2.8158,5,-1,0.3213,-0.3112))
we obtain this results:
summary(b2)
Maximum Likelihood estimation
Newton-Raphson maximisation, 16 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -38.11285
5 free parameters
Estimates:
Estimate Std. error t value Pr( t)
[1,] -2.815780.43548 -6.4660 1.007e-10 ***
[2,] 50.50024 13.17046 3.8344 0.0001259 ***
[3,] -11.533443.31075 -3.4836 0.0004947 ***
[4,] 0.321300.42978 0.7476 0.4547096
[5,] -0.311210.23245 -1.3388 0.1806280
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
NB: a is a binary time series
if we try to estimate the same function using the last version of maxLik
under R.3.0 (64 bit latest version) the estimation do not converge and this
is the error message:
Iteration 15
Parameter:
[1] -2.8146429 51.3042554 -11.7373313 0.3245214 -0.3125767
Gradient:
[1] NaN NaN NaN NaN NaN
Errore in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad,
hessOrig = hess, :
NA in gradient
What causes this?
It could be that the NaNs in the gradients are caused by rounding
errors and/or approximation errors in the numerical
(finite-difference) derivatives when using R.3.0 64 bit, because
different hardware and different software versions (e.g. R,
mathematical libraries, OS) could lead to different rounding errors.
In this case, the specification of a function that returns analytical
gradients could solve the problem.
If this does not solve the problem and you cannot find out the reason
for the NaNs in the analytical gradients yourself, please provide a
reproducible example so that we could help you with this.
Please note that you could also ask questions regarding the maxLik
package via a forum at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
... and please do not forget to cite maxLik in your publications :-)
Best regards,
Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] Stochastic Frontier: Finding the optimal scale/scale efficiency by frontier package
Dear Miao On 25 April 2013 03:26, jpm miao miao...@gmail.com wrote: I am trying to find out the scale efficiency and optimal scale of banks by stochastic frontier analysis given the panel data of bank. I am free to choose any model of stochastic frontier analysis. The only approach I know to work with R is to estimate a translog production function by sfa or other related function in frontier package, and then use the Ray 1998 formula to find the scale efficiency. However, as the textbook Coelli et al 2005 point out that the concavity may not be satisfied, one needs to impose the nonpositive definiteness condition so that the scale efficiency 1. It might be that the true technology is not concave and that the elasticity of scale is larger than one. Indeed, most empirical studies find increasing returns to scale (in many different sectors). Therefore, it is probably inappropriate to impose concavity. How can I do it with frontier package? The frontier package cannot impose concavity on a Translog production function and I am not aware of any software that can do this in a stochastic frontier estimation -- probably, because imposing concavity usually does not make sense. Is there any other SFA model/function in R recommended to find out the scale efficiency and optimal scale? I suggest to plot the elasticity of scale against the firm size. If the elasticity of scale decreases with firm size, then the most productive firm size is at the firm size, where the elasticity of scale is one. However, there are some problems with using the Translog production function (and the Translog distance function) for determining the optimal firm size [1]. [1] http://econpapers.repec.org/RePEc:foi:wpaper:2012_12 If you have further questions regarding the frontier package, I suggest that you use the help forum at frontier's R-Forge site [2]. [2] https://r-forge.r-project.org/projects/frontier/ ... and please do not forget to cite the R packages that you use in your analysis in your publications. Thanks! Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] FRONTIER
Dear Giovanna On 15 February 2013 11:50, Giovanna Ottaviani g...@skogoglandskap.no wrote: Anyone familiar with the package frontier? I have some general questions on how to approach the model design. The R package frontier is based on FRONTIER 4.1 [1] but it has some improvements [2]. The models that can be estimated by the R package frontier are the same as the models that can be estimated by FRONTIER 4.1. These models are explained in the file Front41.pdf that is included in the archive FRONT41-xp1.zip, which is available at [1]. If you have any further questions regarding the frontier package, please use a forum at the package's R-Forge site [3]. [1] http://www.uq.edu.au/economics/cepa/frontier.php [2] http://frontier.r-forge.r-project.org/ [3] http://r-forge.r-project.org/projects/frontier/ Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Error message cs_lu(A) failed: near-singular A (or out of memory)
Dear Rui If you impose the homogeneity (adding-up) restrictions, your system becomes singular, because the error terms of the share equations always sum up to zero. Therefore, you can arbitrarily delete one of the share equations and obtain the coefficients that were not directly estimated by the homogeneity restrictions. Furthermore, you can impose the homogeneity restriction at each single equation by normalization with one input price (numéraire). Finally, I suggest to impose the cross-equation restrictions by the argument restrict.regMat rather than by argument restrict.matrix, because the documentation says the advantage of imposing restrictions on the coefficients by 'restrict.regMat' is that the matrix, which has to be inverted during the estimation, gets smaller by this procedure, while it gets larger if the restrictions are imposed by 'restrict.matrix' and 'restrict.rhs'. I will send you my lecture notes on econometric production analysis with R by private mail. Please do not forget to cite R and the R packages that you use in your publications. If you have further questions regarding system estimation, microeconomic analysis, or stochastic frontier (efficiency) analysis with R, you can use a forum at the R-Forge sites of the systemfit [1,2], micEcon [3,4], or frontier [5] packages/projects, respectively. [1] http://www.systemfit.org/ [2] http://r-forge.r-project.org/projects/systemfit/ [3] http://www.micEcon.org/ [4] http://r-forge.r-project.org/projects/micecon/ [5] http://r-forge.r-project.org/projects/frontier/ Best (holiday) wishes from Copenhagen, Arne On 9 December 2012 23:31, Rui Neiva ruiqne...@gmail.com wrote: Hi there everyone, I have the following model (this is naturally a simplified version just for showing my problem, in case you're wondering this is a translog cost function with the associated cost share equations): C ~ á + â1 log X + â2 log Y + ã1 log Z + ã2 log XX C1 ~ â1 + â2 log YY + ã1 log ZZ Then I have some restrictions on the coefficients, namely that the sum of â equal 1 and the some of ã equal zero So, I've added the following equations to the model C2 ~ 0 + â1.Y1 + â2.Y2 C3 ~ 0 + ã1.Y3 + ã2. Y4 I've created columns in my data frame with values of 0 for variable C3 and values of 1 for Y1, Y2, Y3, Y4 and C2 I'm using the systemfit package to solve a multiple equation system using the SURE method, and using a matrix to impose the restrictions on the coefficients (i.e., that the â1 in all equations is the same value, and the same for all the other coefficients). When I try to run the model without the restricting equations (C2, C3) it runs just fine, but when I add these two equations I get the error: Error in solve(xtOmegaInv %*% xMat2, xtOmegaInv %*% yVec, tol = solvetol) : cs_lu(A) failed: near-singular A (or out of memory) Any ideas on what the problem might be? All the best, Rui Neiva P.S.: I've also posted this question on the Matrix help forum, but since I do not know how active that forum is I've decided to see if anyone in the mailing list would be able to help. [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Coefficient of determination for non-linear equations system (nlsystemfit)
Dear Antti On 20 November 2012 10:24, Antti Simola asimol...@gmail.com wrote: I have estimated system of three linear equations with one non-linear restrictions with nlsystemfit. Please read the FAQ at http://www.systemfit.org/ I was wondering how I can calculate the R-squared (or some alternative coefficient of determination) for the whole system. This is automatically given by linear systemfit but not by nlsystemfit. I can get the values for each of the equations separately, but apparently not for the whole system. You have to do this manually, e.g. by equation (43) in the paper: http://www.jstatsoft.org/v23/i04/paper I'm also wondering why separate equations' objects all appear in together with a command like: nl.system - nlsystemfit() nl.system$eq[i] but e.g. the following produces NULL as the value for the individual objects, e.g. R-squared As the eq component is a list (see documentation), you must use double brackets: nl.system$eq[[i]] eq.1 - nl.system$eq[1] eq.1$r2 You can directly use: nl.system$eq[[i]]$r2 Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] extracting and restricting coefficients
Hi Dereje On 17 October 2012 04:23, Dereje Bacha d_ba...@yahoo.com wrote: I want to fit two equations simultaneously EQ1-Y1~X1+X2 EQ2-Y2~X1+X2 eqsystem-list(Y1HAT=EQ1,Y2HAT=EQ2) fitols-systemfit(eqsystem, method=OLS, data=BB) How do I get coefficients for the first equation? R code coef(fitols$eq[[1]]) How do I restrict coefficient of X2 in the first equation (say , restrict it to less than zero). R code systemfit() can only impose equality constraints. More information about systemfit is available here: http://www.jstatsoft.org/v23/i04/paper and http://www-systemfit.org/ Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] System/SUR equations and panel
Hi Alok On 3 October 2012 17:28, Alok K Bohara, PhD boh...@unm.edu wrote: I was wondering if there is any R package to estimate a two-or three equation system in a panel setting. Say, I have 100 firms and their expenditure on wind generated power and solar-generated power. S_i = X1_i*b1 + u1_i W_i = X2_i* b2 + u2_i (X1 and X2 need not be the same.. but for simplicity, they could be assumed to be similar) Now, I have this firm data set for 10 years... I looked at systemfit package and its investment example and the way it uses plm package was not quite clear. Any help will be much appreciated. Unfortunately, systemfit cannot (yet) apply panel data techniques (e.g. random effects or fixed effects) to estimate a system of equations. However, you can manually estimate a fixed effects model if you add dummy variables for individuals and/or time periods or estimate the model with within transformed data. Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Constraining parameters using tag() in SUR model - ZELIG
On 13 September 2012 18:38, Samantha Azzarello samantha.azzare...@cmegroup.com wrote: Thanks for the help. Ill make sure to cite systemfit along with Zelig. :-) I cannot see how to constrain parameters based on the manual when there is more than 2 eqs using the M Matrix. I have 10 eqns with 6 parameters each and need to constrain all 6 across the 10 equations. For example B1s need to be equal in all 10 eqns, B2 need to be equal in all 10 eqs stc. Can you assist in a brief example? A simple example with three equation and two parameters in each equation: eq1: y1 = a0 + a1 * x1 eq2: y2 = b0 + b1 * x1 eq3: y3 = c0 + c1 * x1 with restrictions a0 = b0 = c0 and a1 = b1 = c1 In fact, you have 4 restrictions (number of equality signs): a0 = b0 b0 = c0 a1 = b1 b1 = c1 Given that you have 4 restrictions and 6 parameters, argument restriction.matrix should be: m - matrix( 0, nrow = 4, ncol = 6 ) # a0 - b0 = 0 m[ 1, 1 ] - 1 m[ 1, 3 ] - -1 # b0 - c0 = 0 m[ 2, 3 ] - 1 m[ 2, 5 ] - -1 # a1 - b1 = 0 m[ 3, 2 ] - 1 m[ 3, 4 ] - -1 # b1 - c1 = 0 m[ 4, 4 ] - 1 m[ 4, 6 ] - -1 Of course, the above restrictions are identical to: a0 = b0 a0 = c0 a1 = b1 a1 = c1 m - matrix( 0, nrow = 4, ncol = 6 ) # a0 - b0 = 0 m[ 1, 1 ] - 1 m[ 1, 3 ] - -1 # a0 - c0 = 0 m[ 2, 1 ] - 1 m[ 2, 5 ] - -1 # a1 - b1 = 0 m[ 3, 2 ] - 1 m[ 3, 4 ] - -1 # a1 - c1 = 0 m[ 4, 2 ] - 1 m[ 4, 6 ] - -1 Alternatively, you can specify the restrictions symbolically using argument restrict.matrix: m - c( eq1_(Intercept) - eq2_(Intercept) = 0, eq2_(Intercept) - eq3_(Intercept) = 0, eq1_x1 - eq2_x1 = 0, eq2_x1 - eq3_x2 = 0 ) Finally, you might use argument restrict.regMat. As you have 6 unconstrained parameters (a0,a1,b0,b1,c0,c1) and two constrained parameters (r0,r1), argument restrict.regMat should be: m - matrix( 0, nrow = 6, ncol = 2 ) m[1,1] - 1 # r0 = a0 m[3,1] - 1 # r0 = b0 m[5,1] - 1 # r0 = c0 m[2,2] - 1 # r1 = a1 m[4,2] - 1 # r1 = b1 m[6,2] - 1 # r1 = c1 All approaches are mathematically equivalent but the latter approach (using argument restrict.regMat) is computationally/numerically preferable. ATTENTION: all of the code above is untested! I am using this in conjunction with a few Perl porgrams already written - but I am open to switching everything I have directly to systemfit if I can get teh parameters constrained OK. /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Constraining parameters using tag() in SUR model - ZELIG
Dear Samantha I suggest that you directly use the systemfit() command rather than the wrapper function zelig(). You can impose parameter restrictions with the arguments restrict.matrix, restrict.rhs, and restrict.regMat. You can use argument restrict.matrix to specify the restrictions numerically (as a matrix) or symbolically (with character strings). The (symbolical and numerical) imposition of parameter restrictions is described in [1], page 18-19 and in the documentation of systemfit() (with examples in the Example section). Please do not forget to cite systemfit in your publications--no matter whether you directly use it or through zelig(). Thanks! [1] http://www.jstatsoft.org/v23/i04/paper /Arne On 12 September 2012 20:05, Samantha Azzarello samantha.azzare...@cmegroup.com wrote: Hello all, I am following some standard code from Zelig manual when using a SUR (Seemingly Unrelated Regression Model) to constrain parameters across equations. Please see code below: setwd(C:/Research/Economics/SUR_FX/Model) # Seemingly Unrelated Regression # Load our library. library(Zelig) library(systemfit) library(zoo) # data. factors - read.table(./RFactors.txt, header = TRUE) returns - read.table(./RReturns.txt, header = TRUE) myData - c(factors,returns) # This is our system of equations. mySys - list(mu1 = USDEUR ~ USDRATE + tag(USDYC, USDYC)+ USDCC + EURRATE + EURYC + EURLY, mu2 = USDGBP ~ USDRATE + tag(USDYC, USDYC) + USDCC + GBPRATE + GBPYC + GBPLY, mu3 = USDCHF ~ USDRATE + tag(USDYC, USDYC) + USDCC + CHFRATE + CHFYC + CHFLY, mu4 = USDSEK ~ USDRATE + tag(USDYC, USDYC) + USDCC + SEKRATE + SEKYC + SEKLY, mu5 = USDNOK ~ USDRATE + tag(USDYC, USDYC) + USDCC + NOKRATE + NOKYC + NOKLY, mu6 = USDJPY ~ USDRATE + tag(USDYC, USDYC) + USDCC + JPYRATE + JPYYC + JPYLY, mu7 = USDSGD ~ USDRATE + tag(USDYC, USDYC) + USDCC + SGDRATE + SGDYC + SGDLY, mu8 = USDAUD ~ USDRATE + tag(USDYC, USDYC) + USDCC + AUDRATE + AUDYC + AUDLY, mu9 = USDCAD ~ USDRATE + tag(USDYC, USDYC) + USDCC + CADRATE + CADYC + CADLY, mu10 = USDNZD ~ USDRATE + tag(USDYC, USDYC) + USDCC + NZDRATE + NZDYC + NZDLY) # Here is our zelig function call. z.out - zelig(mySys,sur, myData) -- After this I am getting error: Error in eval(expr, envir, enclos) : could not find function tag there exists no info (I could find anyways..) RE installing tag() seperate, also tag should be able to be used for the SUR model... Any help would be appreciated. Thanks -- View this message in context: http://r.789695.n4.nabble.com/Constraining-parameters-using-tag-in-SUR-model-ZELIG-tp4642927.html Sent from the R help mailing list archive at Nabble.com. __ R-help@r-project.org mailing list 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] SYSTEMFIT HELP
Dear Arunima On 14 August 2012 08:58, Arunima Haldar arunima.hal...@gmail.com wrote: I want to know whether systemfit can solve simultaneous equations using panel data. systemfit cannot *solve* systems of (simultaneous) equations but it can *estimate* systems of (simultaneous) equations. Do you want to *estimate* a system of (simultaneous) equations? Your question regarding the *estimation* of panel data models is answered in the FAQ section on systemfit's website: http://www.systemfit.org/ If you have further questions regarding the systemfit package, you can either use this mailing list or a forum or tracker on systemfit's R-Forge site: http://r-forge.r-project.org/projects/systemfit/ Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Tobit Fixed Effects
Hi Felipe On 22 March 2012 10:13, Felipe Nunes felipnu...@gmail.com wrote: I'm using censReg() to run a random effects model, but R is having an error. Can you help me understanding why? When I run this model, everything works fine: tob1 - censReg(transfers.cap ~ factor(year) + district2 + gdp.cap - 1, left=0, right=Inf, method=BFGS, nGHQ=15, iterlim=1, data = d3) But this one, does not: tob2 - censReg(transfers.cap ~ factor(year) + I(constituency.coa.v/100) + district2 - 1, left=0, right=Inf, method=BFGS, nGHQ=15, iterlim=1, data = d3) The error is: Error in solve.default(OM) : system is computationally singular: reciprocal condition number = 4.41531e-17 Did you solve this problem in the mean time? If not: I could take a look at it if you send me a reproducible example. Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Systemfit with structural equations and cross equation parameter interaction
Dear Waseem On 7 April 2012 09:20, waseem ahmad waseem...@yahoo.com wrote: I want to estimate simultaneous equation model with panel data. The model looks as follows Y1=a0+a1*X1+a2*X2 Y2=b0+b1*X2+b2*X1 X1=Z1-(Y1/a1) X2=Z2-(Y2/b1) I In this model Y1, Y2, X1 and X2 are endogenous variables; Z1, Z2 are exogenous variables and a0, a1, a2, b0, b1 and b2 are parameters. Could any one please help me how to estimate this model in R. Thanking you in anticipation Please note that your third and fourth equations are non-linear in parameters. For information on estimating systems of non-linear equations with the systemfit package, please take a look at the FAQ section of systemfit's website: http://www.systemfit.org/ Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] How to do 2SLS in R
On 21 March 2012 09:32, Achim Zeileis achim.zeil...@uibk.ac.at wrote: On Wed, 21 Mar 2012, priya.s...@zycus.com wrote: Hi List I want to carry out structural mode. Following Example l have taken from Basic Econometrics- Damodar Gujarati : Look at the systemfit package. For an introduction see the accompanying JSS paper: http://www.jstatsoft.org/v23/i04/ @Achim: Thanks for referring to the systemfit package and the corresponding JSS paper. @Priya: The command tsls (package sem) can only estimate single-equation models, while systemfit can also estimate multiple-equation models. I hope that the documentation of systemfit as well as the JSS paper are sufficiently clear. If you still have questions regarding systemfit, please use a forum or tracker on systemfit's R-Forge site (see http://systemfit.org). Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Demographic Variables in AIDS (Demand System)
Dear xiaodao On 7 March 2012 05:19, xiaodao yanghaow...@utk.edu wrote: I am using aidsEst( ) in micEconAids package to estimate Demand system. But I would like to add more demographic variables in demand system. How can I add those information? for example: mydata-data.frame(p1,p2,p3,p4, s1,s2,s3,s4, totalexp, house,age,urban) priceNames-c(p1,p2,p3,p4) shareNames-c(s1,s2,s3,s4) estResult - aidsEst(priceNames,shareNames,'totalexp',data=mydata, method = LA) for this step, estResult - aidsEst(), how to add more demographic variables in this system? You can use argument shifterNames for this. BTW: You can also ask questions, report problems, etc. about the micEconAids package via a forum or tracker at the R-Forge site [1] of the micEcon project [2]. [1] https://r-forge.r-project.org/projects/micecon/ [2] http://www.micEcon.org/ Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Tobit Fixed Effects
Dear Felipe On 29 September 2011 14:10, Arne Henningsen arne.henning...@googlemail.com wrote: Hi Felipe On 25 September 2011 00:16, Felipe Nunes felipnu...@gmail.com wrote: Hi Arne, my problem persists. I am still using censReg [version - 0.5-7] to run a random effects model in my data (50,000 cases), but I always get the message. tob7 - censReg(transfers.cap ~ pt.pt + psdb.pt + pt.opp + pt.coa + psdb.coa + pib.cap + transfers.cap.lag + pib.cap + ifdm + log(populat) + mayor.vot.per + log(bol.fam.per+0.01) + factor(uf.name) + factor(year) - 1, left=0, right=Inf, method=BHHH, nGHQ=8, iterlim=1, data = pdata2) Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in the initial gradient If I sent you my data set, could you try to help me? I have been struggling with that for two months now. Thanks for sending me your data set. With it, I was able to figure out, where the NAs in the (initial) gradients come from: when calculating the derivatives of the standard normal density function [d dnorm(x) / d x = - dnorm(x) * x], values for x that are larger than approximately 40 (in absolute terms) result in so small values (in absolute terms) that R rounds them to zero. Later, these derivatives are multiplied by some other values and then the logarithm is taken ... and multiplying any number by zero and taking the logarithms gives not a finite number :-( When *densities* of the standard normal distribution become too small, one can use dnorm(x,log=TRUE) and store the logarithm of the small number, which is much larger (in absolute terms) than the density and hence, is not rounded to zero. However, in the case of the *derivative* of the standard normal density function, this is more complicated as log( d dnorm(x) / d x ) = log( dnorm(x) ) + log( - x ) is not defined if x is positive. I will try to solve this problem by case distinction (x0 vs. x0). Or does anybody know a better solution? Finally(!), I have implemented this solution in the censReg() package. Some initial tests (including your model and data) show that the revised calculation of the gradient of the random effects panel data model for censored dependent variables is much more robust to rounding errors. The improved version of the censReg package is not yet via CRAN, but it is available at R-Forge: https://r-forge.r-project.org/R/?group_id=256 If you have further questions or feedback regarding the censReg package, please use a forum or tracker on the R-Forge site of the sampleSelection project: https://r-forge.r-project.org/projects/sampleselection/ Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Prediction from censReg?
Dear Zack On 4 December 2011 03:08, z2.0 zack.abraham...@gmail.com wrote: Wanted to ask if there's an easy way to use 'predict' with objects of class 'censReg', 'maxLik', 'maxim' or 'list'. Have a left-censored dataset, attempting to use a Tobit model and am working with the censReg package. I like how easy it is to move from glm models to predictions with 'predict' and wanted to ask if there was a way to do so with objects output from 'censReg'. Unfortunately, I haven't had time yet to implement a predict() method for censReg objects. Sorry! However, if you have (or anybody else has) time and skills to do this, please do not hesitate to become a project member at R-Forge and implement this method. Thanks! Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Unable to reproduce Stata Heckman sample selection estimates
On 25 November 2011 04:37, Yuan Yuan y.y...@vt.edu wrote: Hello, I am working on reproducing someone's analysis which was done in Stata. The analysis is estimation of a standard Heckman sample selection model (Tobit-2), for which I am using the sampleSelection package and the selection() function. I have a few problems with the estimation: 1) The reported standard error for all estimates is Inf ... vcov(selectionObject) yields Inf in every cell. 2) While the selection equation coefficient estimates are almost exactly the same as the Stata results, the outcome equation coefficient estimates are quite different (different sign in one case, order of magnitude difference in some other cases). 3) I can't seem to figure out how to specify the initial values for the MLE ... whatever argument I pass to start (even of the form coef(selectionObject)), I get the following error: Error in gr[, fixed] - NA : (subscript) logical subscript too long I have to admit I am pretty confused by #1, I feel like I must be doing something wrong, missing something obvious, but I have no idea what. I figure #2 might be because the algorithms (selection and Stata) are just finding different local maxima, but because of #3 I can't test that guess by using different initial values in selection. Let me know if I should provide any more information. Thanks in advance for any pointers in the right direction. Yes, please provide more information (see also the posting guide [1]), e.g. which version of R and which version of the sampleSelection package are you using? Do you estimate the model by the two-step approach or by the 1-step maximum likelihood method? Which commands did use use? Can you send us a reproducible example? Have you read the paper about using the sampleSelection package [2]? [1] http://www.r-project.org/posting-guide.html [2] http://www.jstatsoft.org/v27/i07 Best wishes from copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Estimating model parameters for system of equations
Dear Louise On 15 November 2011 19:03, lstevenson louise.steven...@lifesci.ucsb.edu wrote: Hi all, I'm trying to estimate model parameters in R for a pretty simple system of equations, but I'm having trouble. Here is the system of equations (all derivatives): eqAlgae - (u_Amax * C_A) * (1 - (Q_Amin / Q_A)) eqQuota - (p_max * R_V) / (K_p + R_V) - ((Q_A-Q_Amin)*u_Amax) eqResource - -C_A * (p_max * R_V) / (K_p + R_V) eqSystem - list(C_A = eqAlgae, Q_A = eqQuota, R_V = eqResource) I want to estimate u_Amax, Q_Amin, p_max and Q_Amin with the data I've collected using least squares. I've tried using systemfit but I'm not sure how to write out the equations (my attempt is above but that doesn't work since I haven't given values to the parameters I'm trying to estimate - should I give those parameters initial values?). I've looked into the other functions to get least squares estimates (e.g. lm() ) but I'm not sure how to use that for a system of equations. I have some experience with R but I'm a novice when it comes to parameter estimation, so any help would be much appreciated! Thank you! Your system of equations is non-linear in parameters. As lm() and systemfit() can only estimate models that are linear in parameters, you cannot use these commands to estimate your model. The systemfit package includes the function nlsystemfit() that is intended to estimate systems of non-linear equations. However, nlsystemfit() is still under development and often has convergence problems. Therefore, I wouldn't use it for serious applications. You can estimate your non-linear equations separately with nls(). If you want to estimate your equations jointly, I am afraid that you either have to switch to another software or have to implement the estimation yourself. You could, e.g., minimize the determinant of the residual covariance matrix with optim(), nlm(), nlminb(), or another optimizer or you could maximize the likelihood function of the FIML model using maxLik(). Sorry that I (and R) cannot present you a simple solution! Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
[R] [R-pkgs] mvProbit -- Multivariate Probit Models
Dear R users, I am happy to announce the initial release of the mvProbit package on CRAN (version 0.1-0). This package provides tools for econometric analysis with Multivariate Probit Models. While these models can be estimated also by several other statistical software packages (e.g. LIMDEP/NLOGIT, STATA), mvProbit is much more flexible and powerful in calculating marginal effects. To my best knowledge, mvProbit is the only statistical software package that can calculate various marginal effects including their standard errors at arbitrary user-defined values of the explanatory (and dependent) variables (e.g. at all observations or at the sample mean): marginal effects on the unconditional expectations of the dependent variables and marginal effects on the conditional expectations of each dependent variable at all possible combinations of the other dependent variables. Feedback is very welcome! /Arne -- Arne Henningsen http://www.arne-henningsen.name ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@r-project.org mailing list 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.
Re: [R] sum of functions
On 6 October 2011 16:20, Dimitris.Kapetanakis dimitrios.kapetana...@gmail.com wrote: Dear all, I would like to create a code for semiparametric Klein and Spady's estimator. For that I created a function that provides the log-likelihood function for each observation (so it is a function of betas and i, where i denotes the observation). Now, in order to maximize the log-likelihood function, I have to sum these log-likelihood functions for each i and so to get another function that is a function only of betas and so to maximize it through maxLik for instance. Is that possible? In order to be more clear I give an example of how it could be: Prob1 - function(b, i) g.yj(b,y=1,h.np,i)/(g.yj(b,y=1,h.np,i)+g.yj(b,y=0,h.np,i)) loglik.i- function(b, i) Y[i,]*log(Prob1(b,i))+(1-Y[i,])*log(1-Prob1(b,i)) where b denotes the betas, i the observations, Y is the response vector and g.yj(b,1,h.np,i) a function that I created previously, Prob1(b,i) is a function that gives the conditional probability for observation i and loglik.i(b,i) is function that gives the log-likelihood for observation i. How can I sum the loglik.i(b,i) for each i and remain as a function of b ONLY in order to maximize it??? For exemple this could be done manually by loglik- function(b) loglik.i(b,1)+loglik.i(b,2)+loglik.i(b,3)+….+loglik.i(b,N) but how can I do it automatically for all observations? loglik - function(b) sapply( 1:N, loglik.i, b = b ) Please note that logLik( b ) returns a *vector* of the likelihood contributions of each observation. maxLik() takes the sum of the elements of this vector automatically. If logLik( b ) returns a vector of the likelihood contributions of each observation (rather than just the sum), the BHHH optimisation method can be used. /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Tobit Fixed Effects
Hi Felipe On 25 September 2011 00:16, Felipe Nunes felipnu...@gmail.com wrote: Hi Arne, my problem persists. I am still using censReg [version - 0.5-7] to run a random effects model in my data (50,000 cases), but I always get the message. tob7 - censReg(transfers.cap ~ pt.pt + psdb.pt + pt.opp + pt.coa + psdb.coa + pib.cap + transfers.cap.lag + pib.cap + ifdm + log(populat) + mayor.vot.per + log(bol.fam.per+0.01) + factor(uf.name) + factor(year) - 1, left=0, right=Inf, method=BHHH, nGHQ=8, iterlim=1, data = pdata2) Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in the initial gradient If I sent you my data set, could you try to help me? I have been struggling with that for two months now. Thanks for sending me your data set. With it, I was able to figure out, where the NAs in the (initial) gradients come from: when calculating the derivatives of the standard normal density function [d dnorm(x) / d x = - dnorm(x) * x], values for x that are larger than approximately 40 (in absolute terms) result in so small values (in absolute terms) that R rounds them to zero. Later, these derivatives are multiplied by some other values and then the logarithm is taken ... and multiplying any number by zero and taking the logarithms gives not a finite number :-( When *densities* of the standard normal distribution become too small, one can use dnorm(x,log=TRUE) and store the logarithm of the small number, which is much larger (in absolute terms) than the density and hence, is not rounded to zero. However, in the case of the *derivative* of the standard normal density function, this is more complicated as log( d dnorm(x) / d x ) = log( dnorm(x) ) + log( - x ) is not defined if x is positive. I will try to solve this problem by case distinction (x0 vs. x0). Or does anybody know a better solution? /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] function censReg in panel data setting
On 24 September 2011 23:58, felipnunes felipnu...@gmail.com wrote: Hi guys, did you figure out a way to solve your problem, Igor? I'm still trying to understand what's happening to my estimations. My data set is also big (50,000) and I get the same error: Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in the initial gradient I already tried many different specifications and the error is always the same. I'm fitting the following model: tob7 - censReg(transfers.cap ~ pt.pt + psdb.pt + pt.opp + pt.coa + psdb.coa + pib.cap + transfers.cap.lag + pib.cap + ifdm + log(populat) + mayor.vot.per + log(bol.fam.per+0.01) + factor(uf.name) + factor(year) - 1, left=0, right=Inf, method=BHHH, nGHQ=8, iterlim=1, data = pdata2) Did you upgrade to version 0.5-7 (revision = 1207) from R-Forge? /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Tobit Fixed Effects
On 25 September 2011 00:16, Felipe Nunes felipnu...@gmail.com wrote: Hi Arne, my problem persists. I am still using censReg [version - 0.5-7] to run a random effects model in my data (50,000 cases), but I always get the message. tob7 - censReg(transfers.cap ~ pt.pt + psdb.pt + pt.opp + pt.coa + psdb.coa + pib.cap + transfers.cap.lag + pib.cap + ifdm + log(populat) + mayor.vot.per + log(bol.fam.per+0.01) + factor(uf.name) + factor(year) - 1, left=0, right=Inf, method=BHHH, nGHQ=8, iterlim=1, data = pdata2) Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in the initial gradient If I sent you my data set, could you try to help me? I have been struggling with that for two months now. A reproducible example is always very useful. So please send me your data and the R script that reproduces the error. /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Tobit Fixed Effects
0 0 0 ACACRELANDIA-2008 1 0.0 2007-2008 PSDB/Opo 1 0 0 0 0 psdb.coa psdb.opp time transf.log ACACRELANDIA-2003 0 0 1 -2.236723 ACACRELANDIA-2004 0 0 2 -6.907755 ACACRELANDIA-2005 0 0 3 1.459816 ACACRELANDIA-2006 0 0 4 3.486362 ACACRELANDIA-2007 0 1 1 4.933513 ACACRELANDIA-2008 0 1 2 4.682453 tob6 - censReg(transfers.cap ~ pt.pt + psdb.pt + pt.opp + pt.coa + psdb.coa + transfers.cap.lag + pib.cap + ifdm + log(populat) + mayor.vot.per + bol.fam + factor(uf.name) + factor(year), left=0, right=Inf, method=BHHH, nGHQ=15, data = pdata) Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in the initial gradient Please try to use the latest development version of the censReg package, i.e. version 0.5-7 which is available at R-Forge [1] (see [2], [3], and [4]). [1] https://r-forge.r-project.org/R/?group_id=256 [2] http://tolstoy.newcastle.edu.au/R/e15/help/11/09/7037.html [3] http://tolstoy.newcastle.edu.au/R/e15/help/11/09/7288.html [4] http://tolstoy.newcastle.edu.au/R/e15/help/11/09/7307.html /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Tobit Fixed Effects
On 18 September 2011 04:58, Felipe Nunes felipnu...@gmail.com wrote: Do you guys know how could I increase the time limit for the iterations in censReg? I'm receiving the following message whatever method I use: Newton-Raphson maximisation, 150 iterations Return code 4: Iteration limit exceeded. Log-likelihood: -67680.41 on 42 Df You can use argument iterlim to increase the maximum number of iterations, see documentation of maxNR(). /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Tobit Fixed Effects
Dear Felipe On 15 September 2011 09:58, Felipe Nunes felipnu...@gmail.com wrote: My technical problems were: (1) tob3 - censReg(transfers.cap ~ factor(mayor) + vot.weight + vot.per + transfers.cap.lag + pib.cap + cluster(factor(year)), left=0, data = pool) Using 'censReg' my results never converge. Std. errors are not calculated and the coefficients are weird. Did you estimate the model for the pooled data set (i.e. not considering the panel structure)? Or did you create a panel data frame using pdata.frame(), which is available in package plm, and estimate a random-effects model? What do you exactly mean with never converge? Did you try to increase the maximum number of iterations (argument iterlim)? Did you try to use other optimization methods, e.g. BHHH (argument method)? What is the effect of using cluster() in the model (formula) specification? /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Tobit Fixed Effects
On 15 September 2011 20:09, Felipe Nunes felipnu...@gmail.com wrote: I did both ways (pooled and pdata.frame), but in none I got a result. The coefficients are estimated, but not the std. errors. I'm using BFGS method, but I didn't increase the number of iterations yet. Let me try! ... and I highly recommend to use the BHHH method, particularly for the random effects panel data estimation. /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] function censReg in panel data setting
On 14 September 2011 00:36, Arne Henningsen arne.henning...@googlemail.com wrote: Hi Igors On 13 September 2011 13:27, Igors igors.lahanc...@gmail.com wrote: Any success in finding possible solutions for my problem? Somewhat. The calculation of the log-likelihood values is numerically much more robust/stable now. The log-likelihood contributions of some individuals became minus infinity in your model. This was caused by rounding errors as illustrated in the following simplified example: log( exp( a ) + exp( b ) ) If a and b become smaller than approximately -800, exp( a ) and exp( b ) are rounded to zero and the log of their sum (zero) is minus infinity. I have solved this problem by replacing the above calculation by log( exp( a - c ) + exp( b - c ) ) + c with c = max( a, b ) The source code of the improved censReg package is available on R-Forge [1]; R packages will be available on R-Forge [2] probably within one day. [1] https://r-forge.r-project.org/scm/?group_id=256 [2] https://r-forge.r-project.org/R/?group_id=256 Unfortunately, the calculation of the gradients is still not robust but I expect that I can solve this problem in a similar way as I used to solve the problem with the likelihood function itself. I will continue working on this. I have tried to experiment with size of sample and I get really bad picture. I can't get it work even if sample is ~ 1000 obs. And it is way less than I would like to see working, taking into account my full sample size ~ 540 000 obs. I hope that you have a very fast computer -- or a lot of time for waiting many days or even a few weeks. Now I have also improved the numerical stability of the calculation of the *gradients* of the log-likelihood function. Basically, the problem was that in an equation similar to ( exp(a1) * b1 + exp(a2) * b2 ) / exp(d) a1, a2, and d could have values less than -800 so that exp(a1), exp(a2), and exp(d) became zero and hence, the entire equation became zero divided by zero (NaN) in specific cases. As a1, a2, and d are usually of similar size, I could solve this problem by re-arranging the above equation to exp(a1-d) * b1 + exp(a2-d) * b2 I hope that the numerical stability is sufficiently large now so that you can estimate your large model. Please let me know if it works now. Again, the source code of the improved censReg package is available on R-Forge [1]; R packages will be available on R-Forge [2] probably within one day. Please note that you need at least revision 1207. [1] https://r-forge.r-project.org/scm/?group_id=256 [2] https://r-forge.r-project.org/R/?group_id=256 Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] function censReg in panel data setting
Hi Igors On 13 September 2011 13:27, Igors igors.lahanc...@gmail.com wrote: Any success in finding possible solutions for my problem? Somewhat. The calculation of the log-likelihood values is numerically much more robust/stable now. The log-likelihood contributions of some individuals became minus infinity in your model. This was caused by rounding errors as illustrated in the following simplified example: log( exp( a ) + exp( b ) ) If a and b become smaller than approximately -800, exp( a ) and exp( b ) are rounded to zero and the log of their sum (zero) is minus infinity. I have solved this problem by replacing the above calculation by log( exp( a - c ) + exp( b - c ) ) + c with c = max( a, b ) The source code of the improved censReg package is available on R-Forge [1]; R packages will be available on R-Forge [2] probably within one day. [1] https://r-forge.r-project.org/scm/?group_id=256 [2] https://r-forge.r-project.org/R/?group_id=256 Unfortunately, the calculation of the gradients is still not robust but I expect that I can solve this problem in a similar way as I used to solve the problem with the likelihood function itself. I will continue working on this. I have tried to experiment with size of sample and I get really bad picture. I can't get it work even if sample is ~ 1000 obs. And it is way less than I would like to see working, taking into account my full sample size ~ 540 000 obs. I hope that you have a very fast computer -- or a lot of time for waiting many days or even a few weeks. /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] function censReg in panel data setting
Dear Igors On 9 September 2011 22:00, Igors igors.lahanc...@gmail.com wrote: Thank you for your answer. However it doesn't solve my problem fully. The problem is that I have much bigger data set than I sent to you (it was only a small part : 874 obs.). My full data set is 546718 obs. If I try to use censReg on full data set, then it still gives me the same already mentioned error about Na's in the initial gradient. I have sent you an e-mail with full dataset and the code. Thanks! I would really appreciate if you could check how it works and suggest me how to solve this problem. Now, I get the same error message. The problem is that the log-likelihood function becomes minus infinity. I guess that this problem occurs, because a likelihood contribution of (at least) one individual observation is so small that it is rounded zero and hence, its logarithm is minus infinity. However, I have to take a deeper look into this issue before I can give a definite answer and hopefully find a solution. I have noticed that you use iterlim argumet to specify maximum number of iterations. How this argument could affect possibility of obtaining estimates? Using your small data set, the maximization of the (log) likelihood function in censReg() did not converge within 100 iterations. Therefore, I increased the maximum number of iterations from 100 to 200 -- and then the maximization converged :-( /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] function censReg in panel data setting
On 8 September 2011 09:56, Igors igors.lahanc...@gmail.com wrote: Does censReg expect from panel data to be balanced? No. censReg() can estimate models with unbalanced panel data. The estimation in the code that you sent me indeed does not work but if I remove the user-specified starting values (argument start), it works fine. R summary(UpC) Call: censReg(formula = Power ~ Windspeed, left = -Inf, right = 2000, data = PData_In, nGHQ = 4, method = BHHH, iterlim = 200) Observations: Total Left-censored Uncensored Right-censored 874 0847 27 Coefficients: Estimate Std. error t value Pr( t) (Intercept) -462.26676 19.89517 -23.23 2e-16 *** Windspeed188.387961.72492 109.22 2e-16 *** logSigmaMu 4.540530.03352 135.45 2e-16 *** logSigmaNu 5.087330.02706 188.00 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 BHHH maximisation, 131 iterations Return code 2: successive function values within tolerance limit Log-likelihood: -5538.124 on 4 Df /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] function censReg in panel data setting
Dear Igors On 7 September 2011 10:39, Igors igors.lahanc...@gmail.com wrote: However I am still struggling to obtain model estmates. The same code: UpC - censReg(Power ~ Windspeed, left = -Inf, right = 2000,data=PData_In,method=BHHH,nGHQ = 4,start=c(-691.18,186.79,3.9,3.9)) Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in the initial gradient I have tried to change starting values and regressors in the model several times, but I always get the same mentioned error message. How can I make it work? Is this function maxNRCompute() on the last step of calculation (maximization of the ML)? I had an idea that the error could appear since I have huge sample, but I tried to cut it and it still doesn't work. It is hard to figure out the cause of this error without a reproducible example. Is is possible that you send a reproducible example to me? Could it be that there are NAs in the data or something in the panel data specification is not as censReg() expects it? /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] function censReg in panel data setting
On 6 September 2011 07:51, Igors igors.lahanc...@gmail.com wrote: You chose a suitable solution for the first problem (NA in initial gradient). Unfortunately, the documentation of censReg is not very clear regarding starting values of panel tobit models. Please note that you have to specify 4 starting values: intercept, slope parameter, variance of the individual effects (mu), and variance of the general error term (nu). I have already tried to use 4 parameters, here is what I get: UpC - censReg(Power ~ Windspeed, left = -Inf, right = 2000,data=PData_In,method=BHHH,nGHQ = 4,start=c(-691,189,5,5)) Error in censReg(Power ~ Windspeed, left = -Inf, right = 2000, data = PData_In, : argument 'start' must have length 3 Any ideas how to overcome this one? Could this be an error in package or censReg function? Yes, you are right. This is (was) a bug in the censReg function/package. I have fixed the package. The source code of the new version (0.5-6) is available on R-Forge [1]; R packages will be available on CRAN [2] and R-Forge [3] probably within one or two days. [1] https://r-forge.r-project.org/scm/?group_id=256 [2] http://cran.r-project.org/package=censReg [3] https://r-forge.r-project.org/R/?group_id=256 /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] function censReg in panel data setting
Dear Igors On 5 September 2011 23:58, Igors igors.lahanc...@gmail.com wrote: I have a problem estimating Random Effects model using censReg function. small part of code: UpC - censReg(Power ~ Windspeed, left = -Inf, right = 2000,data=PData_In,method=BHHH,nGHQ = 4) Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in the initial gradient ...then I tried to set starting values myself and here is the error what I got: UpC - censReg(Power ~ Windspeed, left = -Inf, right = 2000,data=PData_In,method=BHHH,nGHQ = 4,start=c(-691,189,5)) Error in names(start) - c(colnames(xMat), logSigmaMu, logSigmaNu) : 'names' attribute [4] must be the same length as the vector [3] How can I solve any of these errors? You chose a suitable solution for the first problem (NA in initial gradient). Unfortunately, the documentation of censReg is not very clear regarding starting values of panel tobit models. Please note that you have to specify 4 starting values: intercept, slope parameter, variance of the individual effects (mu), and variance of the general error term (nu). http://cran.r-project.org/web/packages/censReg/vignettes/censReg.pdf Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Problem running maxLik on R
Dear Subhra On 24 May 2011 02:12, Subhra Bhattacharjee bsub...@gmail.com wrote: I am trying to implement a maximum likelihood estimation using maxLik but keep getting this message Error in is.function(func) : 'func' is missing I am new to R and cannot figure out what's wrong. Any advice will be much appreciated. Please do read the posting guide, the documentation of maxLik, the examples included in this documentation, and the following paper: http://dx.doi.org/10.1007/s00180-010-0217-1 Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] help with the maxBHHH routine
Dear Rohit On 3 May 2011 22:53, Rohit Pandey rohitpandey...@gmail.com wrote: Hello R community, I have been using R's inbuilt maximum likelihood functions, for the different methods (NR, BFGS, etc). I have figured out how to use all of them except the maxBHHH function. This one is different from the others as it requires an observation level gradient. I am using the following syntax: maxBHHH(logLik,grad=nuGradient,finalHessian=BHHH,start=prm,iterlim=2) where logLik is the likelihood function and returns a vector of observation level likelihoods and nuGradient is a function that returns a matrix with each row corresponding to a single observation and the columns corresponding to the gradient values for each parameter (as is mentioned in the online help). however, this gives me the following error: *Error in checkBhhhGrad(g = gr, theta = theta, analytic = (!is.null(attr(f, : the matrix returned by the gradient function (argument 'grad') must have at least as many rows as the number of parameters (10), where each row must correspond to the gradients of the log-likelihood function of an individual (independent) observation: currently, there are (is) 10 parameter(s) but the gradient matrix has only 2 row(s) * It seems it is expecting as many rows as there are parameters. So, I changed my likelihood function so that it would return the transpose of the earlier matrix (hence returning a matrix with rows equaling parameters and columns, observations). However, when I run the function again, I still get an error: *Error in gr[, fixed] - NA : (subscript) logical subscript too long* I have verified that my gradient function, when summed across observations gives the same results as the in built numerical gradient (to the 11th decimal place - after that, they differ since R's function is numerical). I am trying to run a very large estimation (1000's of observations and 821 parameters) and all of the other methods are taking way too much time (days). This method is our last hope and so, any help will be greatly appreciated. Please make yourself familiar with the BHHH algorithm and read the documentation of maxBHHH: it says about argument grad: [...] If the BHHH method is used, ‘grad’ must return a matrix, where rows correspond to the gradient vectors of individual observations and the columns to the individual parameters.[...] More information of the maxLik package is available at: http://dx.doi.org/10.1007/s00180-010-0217-1 Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] look for the package of latent class stochastic frontier
Hi hnlidexi On 10 April 2011 02:47, 李德洗 hnlid...@foxmail.com wrote: Dear all, I want to finished my paper by latent class Stochastic Frontier Analysis , but i can not find the package, is there anyone that may help me As far as I know, there is no R function/package that estimates latent class stochastic frontier models out of the box. Then frontier and Benchmarking packages provide tools for (ordinary) stochastic frontier analysis. Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] maximum likelihood accuracy - comparison with Stata
On 28 March 2011 17:08, Peter Ehlers ehl...@ucalgary.ca wrote: On 2011-03-27 21:37, Alex Olssen wrote: Hi everyone, I am looking to do some manual maximum likelihood estimation in R. I have done a lot of work in Stata and so I have been using output comparisons to get a handle on what is happening. I estimated a simple linear model in R with lm() and also my own maximum likelihood program. I then compared the output with Stata. Two things jumped out at me. Firstly, in Stata my coefficient estimates are identical in both the OLS estimation by -reg- and the maximum likelihood estimation using Stata's ml lf command. In R my coefficient estimates differ slightly between lm() and my own maximum likelihood estimation. Secondly, the estimates for sigma2 are very different between R and Stata. 3.14 in R compared to 1.78 in Stata. I have copied my maximum likelihood program below. It is copied from a great intro to MLE in R by Macro Steenbergen http://artsci.wustl.edu/~jmonogan/computing/r/MLE_in_R.pdf Any comments are welcome. In particular I would like to know why the estimate of sigma2 is so different. I would also like to know about the accuracy of the coefficient estimates. Some comments: 1. Since Kmenta is not in the datasets package, you should mention where you're getting it. (I suppose that for economists, it's obvious, but we're not all economists.) I used the version in the systemfit package. 2. I don't believe your Stata value for sigma2 (by which I assume you mean sigma^2). Are you quoting sigma? 3. I can't reproduce your R value of 3.14 for sigma2. I get 3.725391. 4. (most important) There's a difference between the simple ML estimate (which is biased) and R's unbiased estimate for sigma^2. This dataset has 20 cases, so try multiplying your result by 20/17. 5. As to any difference in coefficients, I would guess that the difference is slight (you don't show what it is); it may be due to the fact that optim() produces approximate solutions, whereas in this simple case, an 'exact' solution is possible (and found by R). 6. In case you weren't aware of it: the stats4 package has an mle() function. ... and there is the maxLik package. http://cran.r-project.org/package=maxLik http://www.springerlink.com/content/973512476428614p/ Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] recommended reading for manual maximum likelihood of a system of equations
Hi Alex On 24 March 2011 04:47, Alex Olssen alex.ols...@gmail.com wrote: I am looking for some recommended reading. I want to read up on the estimation systems of linear equations using maximum likelihood? For a comprehensive review of estimating systems of linear equations, see: http://www.jstatsoft.org/v23/i04/ or Greene, Econometric Analysis, 6th ed., 2008, chapter 10. For ML estimation in R, see: http://www.springerlink.com/content/973512476428614p/ For ML estimations of systems of linear equations, see: Greene, Econometric Analysis, 6th ed., 2008, section 13.6.2. Please note that the iterated SUR estimator should converge to the ML estimator. I have a strongly applied bias, I want to be able to do such estimation myself. Reading with examples would be great. Something which also works through a concentrated maximum likelihood estimation would be even better! I want to use a likelihood ratio test to compare some nested models of linear simultaneous equations. I have estimated the systems using systemfit() and nlsystemfit(). systemfit() can be used with lrtest() which automatically obtains the log-likelihood values. Unfortunately one of my models requires nonlinear coefficient restrictions and lrtest() does not seem to be usable with nlsystemfit(). Yes. Please note: in contrast to systemfit(), the function nlsystemfit() is still under development; it has convergence problems rather often and only a few methods have been implemented yet. I decided it would be a good idea to program the likelihood estimation myself. This would solve the problem and would also be a great learning experience. You are invited to implement the ML estimation in the systemfit package. You could join the developer team at R-Forge. Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] new distribution
Dear Graziella On 22 March 2011 16:08, graziella graziella.bona...@unical.it wrote: I need to introduce new distribution for the error components of a frontier model. For this reason, I use the package Frontier 4.1, The software Frontier 4.1 is not related to R. However, the R package frontier is based on Frontier 4.1. but it includes some bug fixes and improvements and it is much user-friendlier. Therefore, I recommend to use the R package frontier rather than Frontier 4.1. but I need to introduce other distribution function than Normal and Half-Normal (that are those given by default). Unfortunately, the R package frontier (as well as Frontier 4.1) can use only the half-normal and the truncated normal distributions for modelling inefficiency. Implementing further distributions for the inefficiency term is on my to-do list but it has only a mediocre priority. Hence, if you want to use other distributions, you have to implement them yourself either in the (FORTRAN) source code of the frontier package or your could start from scratch, write the likelihood functions (and gradient functions) and use, e.g., the maxLik() package for maximising the likelihood. The source code of the frontier package is hosted on R-Forge. I am sorry that I and the frontier package cannot be more helpful right now. Best wishes from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] R beginner - Error in as.vector(x, mode)
Hi Alex! On 10 March 2011 09:46, Alex Olssen alex.ols...@gmail.com wrote: Hi everyone, I am new to R and keep getting the message Error in as.vector(x, mode) while trying to run nlsystemfit. Below is my exact code. The data is in Stata format because I only recently swapped to R and am trying to compare results from Stata to make sure I know what is going on. I have searched google and read sever R-help articles with this error. They all say the problem is to do with attach(). I haven't used attach() as far as I can tell. However it does seem to be a problem with reading the data. I sometimes get a problem with masking. Finally, when I ran the example in the systemfit documentation for nlsystemfit everything went fine. The major structural difference I can see between the codes is they use data() where as I do not. Any advice would be truly appreciated. library(systemfit) library(foreign) auto74 - read.dta(auto74.dta) Please provide a *reproducible* example. We cannot reproduce your example, because we do not have the file auto74.dta. eqWeight - weight ~ a0 + a1*mpg + a2*price + a3*trunk eqLength - length ~ b0 + b1*mpg + b2*price + b3*trunk model - list(eqWeight, eqLength) start.values - c(a0=0, a1=0, a2=0, a3=0, b0=0, b1=0, b2=0, b3=0) nlsystemfit(OLS, model, start.values, data = auto) Error in as.vector(x, mode) : cannot coerce type 'builtin' to vector of type 'any' I guess that the starting values cause this error. The documentation says: startvals: a list of starting values for the coefficients. Note, it should be a *list*. Please note that the nlsystemfit() function (in contrast to the systemfit() function) is still under development and has convergence problems rather often. Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] R beginner - Error in as.vector(x, mode)
Dear Alex On 10 March 2011 20:31, Alex Olssen alex.ols...@gmail.com wrote: I find it hard to provide a reproducible version of this error. When I use the exact same procedure but get data from data() everything works fine. I.e., I do not think the startvals are the problem - in fact I copied the syntax for the startvals directly from the example in the documentation which does work. The following code is very similar to my code above. The important difference is the way the data is obtained. This code works. library(systemfit) library(foreign) data(mtcars) eqDrat - drat ~ a0 + a1*mpg + a2*qsec eqWt - wt ~ b0 + b1*mpg + b2*qsec model - list(eqDrat, eqWt) start.values - c(a0=0.5, a1=0.5, a2=0, b0=0.5, b1=0.5, b2=0) nlsystemfit(OLS, model, start.values, data = mtcars) I think I could solve my problem if I could use my data through data() instead of through read.table() Does anyone know if this is possible? So, the problem seems to come from your data set. As we do not have your data file, we cannot really help you with this. The command for reading the data (e.g. read.dta(), read.table(), or read.csv()) -- of course -- depends on the file format. /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] VECM - GARCH using maxLik()
On 9 February 2011 16:21, Philipp Grueber philipp.grue...@ebs.edu wrote:
I am still relatively new to R currently working with a large set of
(intraday) financial data (about 3.3m observ.). I wish to estimate a vector
error correction model similar to the following:
r = a_0 + a1*r_(t-1) + a2*r_(t-2) + a3*s_(t-1) + a4*s_(t-2) +
a5*c_(t-1) + e
s = a_0 + a1*s_(t-1) + a2*s_(t-2) + a3*r_(t-1) + a4*r_(t-2) +
a5*c_(t-1) + e
(where c=s-r)
Estimating both equations of the VECM individually using OLS, I discovered
ARCH effects in my errors. Thus, I wish to include a GARCH term. I tried to
work with the fGarch package which, as far as I could see, is pretty
inflexible in terms of model specification and, in particular, unable to
consider exogenous variables (i.e. my error correction term or, at a later
stage, a transition function etc.). I have read about the rgarch package but
have not found a way to specify the model I need.
Thus, a colleague and I tried to write our own code using maxLik() -- but we
ran into some problems.
PROBLEM:
The code below does not work, with the error being: “Return code 100:
Initial value out of range.” We have tried various (more or less realistic)
starting values and manipulated the loglik-function but the error has
persisted. In addition, in our attempts we used a for-loop to estimate sig2.
However, for my dataset, I prefer a solution without such loops because I
assume they would take a lot of time.
QUESTION:
Does anybody know how to fix the code? Any hints what the problem might be?
Maybe someone has already found an alternative way to estimate such a GARCH
model using R?
The following data has exactly the same format as the data I am working
with:
library(maxLik)
library(fGarch)
lag-function(x,k){
c(rep(NA,k),x[1:(length(x)-k)])
}
r- as.vector(garchSim(garchSpec(rseed = 1985), n = 102)[,1]) [3:102]
r_1- lag(r,1)[3:102]
r_2-lag(r,2) [3:102]
s-rnorm(102)[3:102]
s_1- lag(s,1)[3:102]
s_2-lag(s,2) [3:102]
c_1-lag(s-r,1) [3:102]
data-as.matrix(cbind(r,r_1,r_2,s_1,s_2,c_1))
As in my original dataset, I lag every parameter individually -- even though
I know that in many situations, there is an easier way to lag the data. In
my original dataset, this is the result of lagging the time series on a
per-day basis. In other words: My dataset comprises multiple days of data
but I do not want to link the first observation of each day with the last
observation of the previous day. Therefore, it would be helpful if any
possible solution would refer to the above data format.
loglik - function(param) {
res - param[1]
a0 - param[2]
a1 - param[3]
a2 - param[4]
a3 - param[5]
a4 - param[6]
a5 - param[7]
omega - param[8]
alpha - param[9]
beta - param[10]
res - r–a0–a1*r_1–a2*r_2–a3*s_1–a4*s_2-a5*c_1
sig2 - numeric(100)
ll - numeric(100)
sig2[1] - 0
for (i in 2:100) {
sig2[i] - omega + alpha*res[i-1]^2 + beta*sig2[i-1]
ll[i] - -1/2*log(2*pi*sig2[i]) - 1/2*res[i]^2/sig2[i]
}
ll
}
est - maxLik(loglik, start=c(.5,.5,.5,.5,.5,.5,.5,.5,.5,.5))
summary(est)
It is because, your loglik() function returns NAs with the initial values:
loglik(c(.5,.5,.5,.5,.5,.5,.5,.5,.5,.5))
[1] 0.00 -1.281288 -1.129472 -1.897470 -2.009189 -1.603896 -1.422840
[8] -1.295774 -1.372789 -1.386749 -1.180154 -2.296334 -1.497573 -1.304747
[15] -1.728566 -2.013113 -1.550148 -1.331233 -1.619898 -1.474260 -1.257694
[22] -2.639963 -1.551604 -1.729651 -1.439579 -2.336881 -2.140721 -1.789380
[29] -1.589067 -1.366046 -1.237783 -1.122181 -1.047587 -1.426377 -1.202710
[36] -1.300351 -1.298771 -1.217437 -1.096333 -3.968950 -1.889532 -1.914252
[43] -1.647921 -1.792587 -2.206678 -1.812484 -1.637345 -1.847598 -1.793146
[50] -1.492001 -1.289086 -1.210335 -1.824306 -1.580044 -1.350120 -3.977663
[57] -2.213723 -2.022477 -1.704629 -1.599071 -1.789774 -1.421933 -2.004710
[64] -1.645671 -1.544476 -1.371070 -2.947443 -1.845836 -1.840481 -1.664951
[71] -1.564801 -1.327754 -1.463221 -1.230692 -1.099063 -2.188419 -1.395483
[78] -3.801520 -1.934664 -1.641021 -1.508269 -1.670946 -1.762663 -1.496547
[85] -1.836904 -1.855429 -1.850838 -1.883408 -1.808358 -1.487548 -1.893714
[92] -1.794318 -1.552698 -1.367103 -1.223706 -1.152354 -1.145350 -2.550406
[99]NANA
Best wishes from Copenhagen,
Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] VECM - GARCH using maxLik()
On 9 February 2011 18:33, Philipp Grueber philipp.grue...@ebs.edu wrote:
thanks for the quick answer sorry for the mistake. Please find a corrected
version of the code below. Unfortunately, the model still does not work –
due to an error I believed to have overcome: “In log(2 * pi * sig2[i]) :
NaNs produced”
You should make sure that 2 * pi * sig2 is always positive, e.g. by
using a different parametrisation for sig2, e.g.
sig2[i] - exp( omega + alpha*res[i-1]^2 + beta*sig2[i-1] )
... and how can I avoid the for-loop?
for (i in 2:99) {
sig2[i] - omega + alpha*res[i-1]^2 + beta*sig2[i-1]
ll[i] - -1/2*log(2*pi*sig2[i]) - 1/2*res[i]^2/sig2[i]
}
I have no idea for sig2 but you could move ll out of the loop:
ll[2:99] - -1/2*log(2*pi*sig2[2:99]) - 1/2*res[2:99]^2/sig2[2:99]
/Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] fast optimization routines in R
Hi Pia! On 7 February 2011 21:44, ppinger pin...@uni-mannheim.de wrote: I am looking for a maximization routine that I can use to maximize a large variety of relatively complex likelihoods. I undertand (from previous posts) that coding the objective function more efficiently can help. However, the optimization routine employed seems important too. So far, I have tried the optimization routines optim, maxlik, trust and nlminb. The latter two are much faster than the first ones but nevertheless, it seems to me as if these routines were rather slow, when compared to some of the optimizers in MATLAB. maxLik() includes several different optimisation routines. Did you try all of them? /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Hausman Test
Hi Holger! On 16 January 2011 15:53, Holger Steinmetz holger.steinm...@web.de wrote: One follow up question. The Hausman-test always gives me a p-value of 1 - no matter how small the statistic is. I now generated orthogonal regressors (X1-X3) and the test gives me Hausman specification test for consistency of the 3SLS estimation data: data Hausman = -0.0138, df = 2, p-value = 1 What is confusing to me is the 3SLS. I am just beginning to learn about instrumental variables (I am a psychologist ;) Perhaps that's a problem? As a background, here's the complete simulation: W = rnorm(1000) X2 = rnorm(1000) X3 = rnorm(1000) X1 = .5*W + rnorm(1000) Y = .4*X1 + .5*X2 + .6*X3 + rnorm(1000) data = as.data.frame(cbind(X1,X2,X3,Y,W)) fit2sls - systemfit(Y~X1,data=data,method=2SLS,inst=~W) fitOLS - systemfit(Y~X1,data=data,method=OLS) print(hausman.systemfit(fitOLS, fit2sls)) Please do read the documentation of hausman.systemfit(). I regret that comparing 2SLS with OLS results has not been implemented yet: == part of documentation of hausman.systemfit() = Usage: hausman.systemfit( results2sls, results3sls ) Arguments: results2sls : result of a _2SLS_ (limited information) estimation returned by ‘systemfit’. results3sls : result of a _3SLS_ (full information) estimation returned by ‘systemfit’. Details: The null hypotheses of the test is that all exogenous variables are uncorrelated with all disturbance terms. Under this hypothesis both the 2SLS and the 3SLS estimator are consistent but only the 3SLS estimator is (asymptotically) efficient. Under the alternative hypothesis the 2SLS estimator is consistent but the 3SLS estimator is inconsistent. The Hausman test statistic is m = ( b_2 - b_3 )' ( V_2 - V_3 ) ( b_2 - b_3 ) where $b_2$ and $V_2$ are the estimated coefficients and their variance covariance matrix of a _2SLS_ estimation and $b_3$ and $V_3$ are the estimated coefficients and their variance covariance matrix of a _3SLS_ estimation. = Please don't hesitate to write a new version of hausman.systemfit() that can also compare 2SLS with OLS results. Best regards from Copenhagen, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Hausman Test
Hi Achim! On 16 January 2011 16:37, Achim Zeileis achim.zeil...@uibk.ac.at wrote: Arne: Unless I'm missing something, hausman.systemfit() essentially does the right thing and computes the right statistic and p-value (see my other mail to Holger). Maybe some preliminary check on the input objects could be used for determining the right order of models. Thanks for the response and the suggestions. Adding a check for the input objects and extending the documentation is a good idea! I will change the systemfit package accordingly in the future. /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Tobit model on unbalanced panel
On 23 November 2010 23:32, Achim Zeileis achim.zeil...@uibk.ac.at wrote: On Tue, 23 Nov 2010, Liang Peng wrote: Appreciate any suggestions regarding how to fit an unbalanced panel data to a Tobit model using R functions. I am trying to analyze how real estate capital expenditures (CapEx) are affected by market conditions using a panel Tobit model. The CapEx is either positive or 0, so it is censored. The data are unbalanced panel, including the CapEx of about 5000 properties over about 40 quarters, with the starting and ending quarters differ across properties. In case you are not familiar with the term Tobit, the model is essentially the following. The true value of the CapEx of property i in quarter t, Y*[i,t], is determined by a property fixed effect (dummy) a[i], explanatory variables x[I,t], and an error: Y*[i,t]=a[i]+b*X[i,t]+u[i,t]. Further, the observed CapEx, Y[i,t], equals Y*[i,t] if Y[i,t]0 and equals 0 otherwise. Now I observe Y[i,t] and X[i,t], and the purpose of this analysis is to estimate coefficients b. Thanks for any suggestions! I think that the censReg package might be worth exploring, see http://CRAN.R-project.org/package=censReg I haven't looked at the package in detail yet, but it says that it is concerned with estimation of censored regression (Tobit) models with cross-section and panel data. Yes, the censReg function/package can estimate random effects tobit models with balanced and unbalanced panel data. However, the estimation with *unbalanced* panel data has not been thoroughly tested yet. If you experience any problems, please report them via the tracker of the sampleSelection project on R-Forge [1]. Thanks! [1] https://r-forge.r-project.org/tracker/?group_id=256 Best regards, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Tobit model on unbalanced panel
Hi Terry On 24 November 2010 15:34, Terry Therneau thern...@mayo.edu wrote: For tobit regression, also see the last example of help(survreg) in the survival package. How does the survreg() function of the survival package account for the panel structure of the data? Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Hausman test for endogeneity
Hi Holger On 10 October 2010 15:36, Holger Steinmetz holger.steinm...@web.de wrote: After inspecting the options, I *guess* that systemfit is what I need. However, I absolutely don't understand how it works. I searched long for a detailed documentation (beyond the rather cryptic standard documentation) but found none. Has anybody references/advises how to conduct the test? A paper describing the systemfit package has been published in the journal of statistical software: http://www.jstatsoft.org/v23/i04/paper It describes the Hausman test for testing the consistency of the 3SLS estimates against the 2SLS estimates (see sections 2.8 and 4.6). I guess (but I am not sure -- maybe others can comment on this) that you test for the endogeneity of regressors, e.g., by fitSur - systemfit( myFormula, data = myData, method = SUR ) fit3sls - systemfit( myFormula, data = myData, method = 3SLS, inst = myInst ) hausman.systemfit( fit3sls, fitSur ) If some regressors are endogenous, the SUR estimates are inconsistent but the 3SLS estimates are consistent given that the instrumental variables are exogenous. However, if all regressors are exogenous, both estimates should be consistent but the SUR estimates should be more efficient. Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Can I monitor the iterative/convergence process while using Optim or MaxLik?
On 14 September 2010 21:49, Sally Luo shali...@gmail.com wrote: Is it possible that I have the estimates from each step/iteration shown on the computer screen in order to monitor the process while I am using Optim or MaxLik? You can call maxLik() with argument print.level set to a value larger than 0. /Arne -- Arne Henningsen http://www.arne-henningsen.name [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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.
Re: [R] How to use lm() output for systemfit() 'Seemingly unrelated regression'
On 3 September 2010 12:03, zbynek.jano...@gmail.com zbynek.jano...@centrum.cz wrote: I am having problem using output of lm() function for further analysing using systemfit package. Basicaly, the problem s following - I generate several formulas using lm() fo1 - lm(r98[,2] ~ f98[,1] + f98[,2] + ... + f98[,43]) fo2 - lm(r98[,1] ~ f98[,1] + f98[,2] + ... + f98[,43]) and than I want to estimate a general model using package systemfit. fitsur - systemfit(SUR,list(as.formula(fo),as.formula(foo))) and I get following error: Error in systemfit(SUR, list(as.formula(fo), as.formula(foo))) : argument 'formula' must be an object of class 'formula' or a list of objects of class 'formula' and I am not able to find where the problem is. The problem remains even when I try: fo1 - r98[,2] ~ f98[,1] + f98[,2] + ... + f98[,43] fo2 - r98[,1] ~ f98[,1] + f98[,2] + ... + f98[,43] instead of using lm() Could someone give me a hand? I am quite new to R, so possibly the solutions is simple:) Yes, the formula(s) must be the *first* argument of the systemfit() command, e.g.: R fitsur - systemfit( list(as.formula(fo),as.formula(foo)), method = SUR ) /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] maxNR in maxLik package never stops
Hi Martin! Which version of maxLik do you use? (You can figure this out, e.g. by typing 'help( package = maxLik )' at the R prompt.) If you are not using version 0.8-0, please upgrade to the latest version. You can get more information on the optimization process by adding print.level=4 as argument of maxLik() or maxNR(). You could send us also an example which allows us to replicate your problem. Best wishes, Arne On 29 August 2010 01:06, Martin Boďa h...@azet.sk wrote: Greetings, I use maxNR function under maxLik package to find the REML estimates of the parameters of variance components in heteroskedastic linear regression models. I assume that in the model there is additive/multiplicative/mixed heteroskedasticity and I need estimate the respective parameters of additive/multiplicative/mixed variance components. For my research purposes I make a batch of simulations and estimate for the given design matrix but different (simulated) response vectors the variance component parameters. However, it happens very frequently, usually after the first simulation that the log-likelihood function maximization step (for which maxNR is used) R freezes. It appears nothing wrong with R, it is working but the maximization never stops (frankly, I cannot say that it is really never, I did not allow it more time than one day). Therefore, my simulations crash - I cannot even perform one simulation. I see two possible solutions: 1. Using some additional parameters of maxNR (such as number of iterations or something). 2. Using some function that could help me control the length of the maximization step - if it were to take to long, it would be evaluated as an unsuccessful simulation. Since I do not have much experience with R, it is clear to me that I cannot resolve the thing myself. Maybe there are other viable solutions. I shall be much obliged for your help and advice. Martin Boďa (Matej Bel University / Faculty of Natural Sciences / Banská Bystrica / Slovakia). [[alternative HTML version deleted]] __ R-help@r-project.org mailing list 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Robust SE Heteroskedasticity-consistent estimation
On 13 May 2010 00:41, RATIARISON Eric eratiari...@monceauassurances.com wrote:
Hi, here my new version:
I submit you my case:
( Pseudo likehood for exponential family with offset )
loglik - function(param) {
b-param
m=as.vector(of+z%*%b)
ll - sum(v*m-exp(m)) }
gradlik - function(param) {
b-param
m=as.vector(of+z%*%b)
gg-(v-exp(m))*z }
hesslik - function(param) {
b-param
m=as.vector(of+z%*%b)
hh - -t(exp(m)*z)%*%z }
resMaxlik -
maxLik(loglik,grad=gradlik,hess=hesslik,start=sc,method=nr,tol=1e-4);
resMaxlik$offset-of
resMaxlik$x-z
resMaxlik$y-v
estfun - function(obj,...)
{
m=as.vector(obj$offset+obj$x%*%obj$estimate)
(obj$y-exp(m))*obj$x
}
M - function(obj, adjust = FALSE, ...)
{
psi - estfun(obj)
k - NCOL(psi)
n - NROW(psi)
rval - crossprod(as.matrix(psi))/n
if(adjust) rval - n/(n - k) * rval
rval
}
B - function(obj, ...)
{ as.matrix(ginv(obj$hessian))
}
So S -B(resMaxlik)*M(resMaxlik)*B(resMaxlik) directly. It's ok.
But I call sandwich function like this : sandwich(resMaxlik,meat.=M,bread.=B)
It returns a error message : Erreur dans UseMethod(estfun) :
pas de méthode pour 'estfun' applicable pour un objet de classe c('maxLik',
'maxim', 'list')
I have added an estfun() method and a bread() method for objects of
class maxLik so that you can use sandwich now.
Please note:
a) The methods can be applied only if maxLik() was called with
argument grad equal to a gradient function or (if no gradient
function is specified) argument logLik equal to a log-likelihood
function that return the gradients or log-likelihood values,
respectively, for *each observation* (i.e. in the BHHH form).
b) The package with the new features is not yet available on CRAN but
on R-Forge: http://r-forge.r-project.org/scm/?group_id=254
Please let me know if you experience any problems.
/Arne
-Message d'origine-
De : Arne Henningsen [mailto:arne.henning...@googlemail.com]
Envoyé : mardi 11 mai 2010 08:25
À : Achim Zeileis; RATIARISON Eric; r-help@r-project.org; Ott-Siim Toomet
Objet : Re: [R] Robust SE Heteroskedasticity-consistent estimation
On 11 May 2010 00:52, Achim Zeileis achim.zeil...@uibk.ac.at wrote:
On Mon, 10 May 2010, RATIARISON Eric wrote:
I'm using maxlik with functions specified (L, his gradient hessian).
Now I would like determine some robust standard errors of my estimators.
So I 'm try to use vcovHC, or hccm or robcov for example
but in use one of them with my result of maxlik, I've a the following
error message :
Erreur dans terms.default(object) : no terms component
Is there some attributes to give to maxlik objet for fitting the call
of vcovHC?
This is discussed in
vignette(sandwich-OOP, package = sandwich)
one of the vignettes accompanying the sandwich package that provides the
vcovHC() function. At the very least, you need an estfun() method which
extracts the gradient contributions per observation. Then you need a bread()
function, typically based on the observed Hessian. Then you can compute the
basic sandwich() estimators.
Is it possible to implement the estfun() method and the bread()
function in the maxLik package so that vcovHC() can be easily used by
all users of the maxLik package? If yes: should we (Eric, Achim, Ott,
Arne) implement this feature together?
/Arne
--
Arne Henningsen
http://www.arne-henningsen.name
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] Robust SE Heteroskedasticity-consistent estimation
On 17 May 2010 22:55, RATIARISON Eric eratiari...@monceauassurances.com wrote:
It's ok Arne, i've build the MS windows binary.
And the result is ok (sandwich function runs)with this next modifications in
my user specified functions
But with the following functions (individual one):
sc= as.vector( c(log(mean(v)),rep(0,dim(z)[2]-1)))
loglik - function(param) {
b-param
m=as.vector(of+z%*%b)
ll - (v*m-exp(m)) } # no sum
gradlik - function(param) {
b-param
m=as.vector(of+z%*%b)
gg -(v-exp(m))*z } #here I've replaced %*% by *
hesslik - function(param) {
b-param
m=as.vector(of+z%*%b)
hh - -t(exp(m)*z)%*%z }
resMaxlik - maxLik(loglik,grad=gradlik,hess=hesslik,start=sc,method=nr)
however, I've a question again about bread function. Does it use the result
of hessian defined by user?
Yes.
When I do all.equal(vcov(resMaxlik),inv(resMaxlik$hessian))
I've this message : [1] Mean relative difference: 2 what does it means?
The covariance matrix is the inverse of the *negative* Hessian. So
all.equal(vcov(resMaxlik),solve(-hessian(resMaxlik)))
should return TRUE
Lastly, in the last official manual Maxlik.pdf, the sentence However, no
attempt is made to correct the resulting variance-covariance matrix
Do you mean this paragraph?
‘maxLik’ support constrained optimization in the sense that
constraints are passed further to the underlying optimization
routines, and suitable default method is selected. However, no
attempt is made to correct the resulting variance-covariance
matrix. Hence, the inference may be wrong. A corresponding
warning is issued by the summary method.
This means that the constraints are not taken into account when maxLik
calculates the covariance matrix of the estimates -- the constraints
are simply ignored.
Should be corrected with possible use of sandwich should'nt it?
I don't think that just using the sandwich method is sufficient to
calculate consistent covariance matrices, when there are (equality or
inequality) constraints on the parameters.
/Arne
-Message d'origine-
De : Arne Henningsen [mailto:arne.henning...@googlemail.com]
Envoyé : lundi 17 mai 2010 17:58
À : RATIARISON Eric
Objet : Re: [R] Robust SE Heteroskedasticity-consistent estimation
On 17 May 2010 17:35, RATIARISON Eric eratiari...@monceauassurances.com
wrote:
Thank you Arne,
I'm working on PC computer at my work with firewall.
How can I reach the new version for testing (in comparing with the gauss
procedure named CML) ?
I don't understand the link joined.
I have attached the package. You can easily install it if you (also) have a
Unix-like OS. If you have MS-Windows, I could build a MS-Windows binary
package on http://win-builder.r-project.org/ tomorrow (or you could change
the maintainer field and do it yourself).
/Arne
-Message d'origine-
De : Arne Henningsen [mailto:arne.henning...@googlemail.com]
Envoyé : lundi 17 mai 2010 17:12
À : RATIARISON Eric
Cc : Achim Zeileis; r-help@r-project.org; Ott-Siim Toomet Objet : Re:
[R] Robust SE Heteroskedasticity-consistent estimation
On 13 May 2010 00:41, RATIARISON Eric eratiari...@monceauassurances.com
wrote:
Hi, here my new version:
I submit you my case:
( Pseudo likehood for exponential family with offset )
loglik - function(param) {
b-param
m=as.vector(of+z%*%b)
ll - sum(v*m-exp(m)) }
gradlik - function(param) {
b-param
m=as.vector(of+z%*%b)
gg-(v-exp(m))*z }
hesslik - function(param) {
b-param
m=as.vector(of+z%*%b)
hh - -t(exp(m)*z)*z }
resMaxlik -
maxLik(loglik,grad=gradlik,hess=hesslik,start=sc,method=nr,tol=1e-4
);
resMaxlik$offset-of
resMaxlik$x-z
resMaxlik$y-v
estfun - function(obj,...)
{
m=as.vector(obj$offset+obj$x%*%obj$estimate)
(obj$y-exp(m))*obj$x
}
M - function(obj, adjust = FALSE, ...) { psi - estfun(obj) k -
NCOL(psi) n - NROW(psi) rval - crossprod(as.matrix(psi))/n
if(adjust) rval - n/(n - k) * rval
rval
}
B - function(obj, ...)
{ as.matrix(ginv(obj$hessian))
}
So S -B(resMaxlik)*M(resMaxlik)*B(resMaxlik) directly. It's ok.
But I call sandwich function like this :
sandwich(resMaxlik,meat.=M,bread.=B)
It returns a error message : Erreur dans UseMethod(estfun) :
pas de méthode pour 'estfun' applicable pour un objet de classe
c('maxLik', 'maxim', 'list')
I have added an estfun() method and a bread() method for objects of
class maxLik so that you can use sandwich now.
Please note:
a) The methods can be applied only if maxLik() was called with
argument grad equal to a gradient function or (if no gradient
function is specified) argument logLik equal to a log-likelihood
function that return the gradients or log-likelihood values,
respectively, for *each observation* (i.e. in the BHHH form).
b) The package with the new features is not yet available on CRAN but
on R-Forge: http://r-forge.r-project.org
Re: [R] Robust SE Heteroskedasticity-consistent estimation
On 11 May 2010 00:52, Achim Zeileis achim.zeil...@uibk.ac.at wrote: On Mon, 10 May 2010, RATIARISON Eric wrote: I'm using maxlik with functions specified (L, his gradient hessian). Now I would like determine some robust standard errors of my estimators. So I 'm try to use vcovHC, or hccm or robcov for example but in use one of them with my result of maxlik, I've a the following error message : Erreur dans terms.default(object) : no terms component Is there some attributes to give to maxlik objet for fitting the call of vcovHC? This is discussed in vignette(sandwich-OOP, package = sandwich) one of the vignettes accompanying the sandwich package that provides the vcovHC() function. At the very least, you need an estfun() method which extracts the gradient contributions per observation. Then you need a bread() function, typically based on the observed Hessian. Then you can compute the basic sandwich() estimators. Is it possible to implement the estfun() method and the bread() function in the maxLik package so that vcovHC() can be easily used by all users of the maxLik package? If yes: should we (Eric, Achim, Ott, Arne) implement this feature together? /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
[R] [R-pkgs] micEcon split into micEconSNQP, micEconCES, and micEcon
In December, the packages miscTools and micEconAids were separated from the micEcon package. Now, two further packages have been separated from the micEcon package: micEconSNQP and micEconCES. The micEconSNQP package (version 0.6-2) provides tools for production analysis with the Symmetric Normalized Quadratic (SNQ) profit function (also known as the Generalized McFadden profit function), e.g. snqProfitEst() for estimating the SNQ profit function, snqProfitCalc() for calculating the endogenous variable of the SNQ profit function, or snqProfitEla() for calculating price elasticities of the SNQ profit function. The micEconCES package (version 0.6-6) provides tools for estimating the Constant Elasticity of Scale (CES) function. Function cesEst() has been tremendously improved. It can estimate the CES function using much more estimation methods (optimizers) now, e.g. the Levenberg-Marquardt or the Differential Evolution algorithm. Furthermore a vignette (supplementary documentation) Estimating the CES Function in R: Package micEconCES (written by Geraldine Henningsen and Arne Henningsen) [1] has been added. The micEcon package (version 0.6-2) includes the remaining parts of the old micEcon package, e.g. for economic analysis with Cobb-Douglas functions, translog functions, quadratic functions, or linearly homogeneous non-parametric functions. All these packages are available on CRAN [2,3,4] and R-Forge [5]. [1] http://cran.r-project.org/web/packages/micEconCES/vignettes/CES.pdf [2] http://cran.r-project.org/package=micEcon [3] http://cran.r-project.org/package=micEconSNQP [4] http://cran.r-project.org/package=micEconCES [5] http://r-forge.r-project.org/projects/micecon/ Feedback is highly welcome! /Arne -- Arne Henningsen http://www.arne-henningsen.name ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@r-project.org mailing list 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.
Re: [R] maxNR - Error in p(a, b) : element 1 is empty; the part of the args list of '*' being evaluated was: (b, t)
On 22 March 2010 16:31, 4-real danielkjae...@gmail.com wrote:
Hello everyone...
We were trying to implement the Newton-Raphson method in R, and estimate the
parameters a and b, of a function, F, however we can't seem to implement
this the right way. Hope you can show me the right way to do this. I think
what we want R to do is to read the data from the website and then peform
maxNR on the function, F. Btw the version of R being used is RGui for
Windows if it helps to know this.
R-code below:
library(maxLik)
require(maxLik)
x -
read.table('http://www.math.ku.dk/kurser/2008-09/blok4/stat2/doku/data/Eksempel_6_3.txt',
header = TRUE);
t - log(x$Koncentration);
X - x$Status;
p - function(a,b) exp(a+b*t)/(1+exp(a+b*t));
S - sum(X);
SP - sum(t*X);
F - function(a,b) {
+ c(sum(p(a,b)) - S,
+ sum(t*p(a,b)) - SP)
+ }
z - maxNR(F, start=1, print.level=2)
Error in p(a, b) : element 1 is empty;
the part of the args list of '*' being evaluated was:
(b, t)
You forgot to provide argument b for function F:
R z - maxNR(F, start=1, b=1, print.level=2)
- Initial parameters: -
fcn value: -15.42843
parameter initial gradient free
[1,] 1-22.490161
Condition number of the (active) hessian: 1
-Iteration 1 -
-Iteration 2 -
-Iteration 3 -
-Iteration 4 -
-Iteration 5 -
-Iteration 6 -
--
gradient close to zero. May be a solution
6 iterations
estimate: -0.6288598
Function value: -1.536789
/Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] Constrained non linear regression using ML
On 17 March 2010 14:22, Gabor Grothendieck ggrothendi...@gmail.com wrote: Contact the maintainer regarding problems with the package. Not sure if this is acceptable but if you get it to run you could consider just dropping the variables from your model that correspond to active constraints. Also try the maxLik package. You will have to define the likelihood yourself but it does support constraints. Yes. And specifying the likelihood function is probably (depending on your distributional assumptions) not too complicated. BTW: Even if your y follows a beta distribution, it does not mean that your error term also follows a beta distribution. And it the distribution of the error term which is crucial for specifying the likelihood function. /Arne On Wed, Mar 17, 2010 at 9:07 AM, Corrado ct...@york.ac.uk wrote: Dear Gabor, 1) The constraints are active, at least from a formal point view. 3) I have tried several times to run betareg.fit on the data, and the only thing I can obtain is the very strange error: Error in dimnames(x) - dn : length of 'dimnames' [2] not equal to array extent The error is strange because, because the function dimnames is not called anywhere. Regards Gabor Grothendieck wrote: Try it anyways -- maybe none of your constraints are active. On Wed, Mar 17, 2010 at 6:01 AM, Corrado ct...@york.ac.uk wrote: Dear Gabor, dear R users, I had already read the betareg documentation. As far as I can understand from the help, it does not allow for constrained regression. Regards Gabor Grothendieck wrote: Check out the betareg package. On Tue, Mar 16, 2010 at 2:58 PM, Corrado ct...@york.ac.uk wrote: Dear R users, I have to fit the non linear regression: y~1-exp(-(k0+k1*p1+k2*p2+ +kn*pn)) where ki=0 for each i in [1 n] and pi are on R+. I am using, at the moment, nls, but I would rather use a Maximum Likelhood based algorithm. The error is not necessarily normally distributed. y is approximately beta distributed, and the volume of data is medium to large (the y,pi may have ~ 40,000 elements). I have studied the packages in the task views Optimisation and Robust Statistical Methods, but I did look like what I was looking for was there. Maybe I am wrong. The nearest thing was nlrob, but even that does not allow for constraints, as far as I can understand. Any suggestion? Regards -- Corrado Topi PhD Researcher Global Climate Change and Biodiversity Area 18,Department of Biology University of York, York, YO10 5YW, UK Phone: + 44 (0) 1904 328645, E-mail: ct...@york.ac.uk __ R-help@r-project.org mailing list 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. -- Corrado Topi PhD Researcher Global Climate Change and Biodiversity Area 18,Department of Biology University of York, York, YO10 5YW, UK Phone: + 44 (0) 1904 328645, E-mail: ct...@york.ac.uk -- Corrado Topi PhD Researcher Global Climate Change and Biodiversity Area 18,Department of Biology University of York, York, YO10 5YW, UK Phone: + 44 (0) 1904 328645, E-mail: ct...@york.ac.uk __ R-help@r-project.org mailing list 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. -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Question regarding to maxNR
Hi Kate!
On Fri, Mar 12, 2010 at 6:20 AM, kate yh...@illinois.edu wrote:
Hi R-users,
Recently, I use maxNR function to find maximizer. I have error appears as
follows
Error in maxNRCompute(fn = fn, grad = grad, hess = hess, start = start, :
NA in the initial gradient
My code is
mu=2
s=1
n=300
library(maxLik)
set.seed(1004)
x-rcauchy(n,mu,s)
loglik-function(mu)
{
log(prod(dcauchy(x,mu,s)))
}
maxNR(loglik,start=median(x))$estimate
Does anyone know how to solve this problem?
Yes :-)
Algebraically, log(prod(z)) is equal to sum(log(z)) but these two
expressions might return different numbers on digital computers:
R ll-function(mu) { sum(log(dcauchy(x,mu,s))) }
R ll(2)
[1] -754.4928
R loglik(2)
[1] -Inf
R maxNR(ll,start=median(x))$estimate
[1] 2.075059
Best wishes from Copenhagen,
Arne
--
Arne Henningsen
http://www.arne-henningsen.name
__
R-help@r-project.org mailing list
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.
Re: [R] Tutorials and scripts of Stochastic Frontier Analysis and Linear Programming.
Hi Marcus! On Sun, Feb 21, 2010 at 1:47 PM, Marcus Vinicius mvi...@gmail.com wrote: I want to program my own models about Stochastic Frontier Analysis and Linear programming (Data Envelopment Analysis). In this context, is there anyone that may help me with some simple tutorials and scripts about these issues? What do you mean with own models? Have you developed new theoretical models that are not (yet) available in current software package? Various types of Stochastic Frontier Analysis (SFA) can be performed by the frontier package [1,2] and there are (at least) three R packages for Data Envelopment Analysis: DEA [3], FEAR [4] (academic use without payments but NOT open source :-( ), and one of my colleagues has an unpublished package for DEA (I can ask him to send it to you). [1] http://cran.r-project.org/package=frontier [2] http://r-forge.r-project.org/projects/frontier/ [3] http://cran.r-project.org/package=DEA [4] http://www.clemson.edu/economics/faculty/wilson/Software/FEAR/fear.html Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Functions for QUAIDS and nonlinear SUR?
On Sat, Jan 9, 2010 at 1:21 AM, Werner W. pensterfuz...@yahoo.de wrote: I would like to estimate a quadratic almost ideal demand system in R which is estimated usually by nonlinear seemingly unrelated regression. But there is no such function in R yet The systemfit package has the function nlsystemfit() for estimating systems of non-linear equations, e.g. by non-linear SUR. However, in contrast to the systemfit() function for estimating systems of linear equations, the function nlsystemfit() is still under development and has convergence problems rather often. So, I cannot recommend using nlsystemfit() for an important analysis. :-( but it is readily available in STATA (nlsur), see B. Poi (2008): Demand-system estimation: Update, Stata Journal 8(4). Now I am thinking, what is quicker learning to program STATA which seems not really comfortable for programming or implement the method in R which might be above my head in terms of econometrics. You do not have to start from scratch but you could improve the nlsystemfit() function, e.g. by implementing analytical gradients of the objective function -- and I could assist you with this. If you are interested in improving nlsystemfit(), please apply at R-Forge [1] for getting write access to systemfit's SVN repository. [1] http://r-forge.r-project.org/projects/systemfit/ /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Interpreting coefficient in selection and outcome Heckman models in sampleSelection
Hi Mark! On Sun, Jan 3, 2010 at 9:08 PM, Mark Bulling mark.bull...@googlemail.com wrote: Hi there Within sampleSelection, I'm trying to calculate the marginal effects for variables that are present in both the selection and outcome models. For example, age might have a positive effect on probability of selection, but then a negative effect on the outcome variable. i.e. Model-selection(participation~age, frequency~age, ...) Documentation elsewhere describes one method for doing this in Stata based on Sigelman and Zeng: http://polisci.osu.edu/prl/Selection%20Models.pdf - see page 16. I'd like to replicate this in r, but wanted to check I'm not reinventing the wheel, before doing so. I don't know a function/method that does this in R. So if you want to implement this in R, I suggest that you add a marginalEffects (or similar) method for objects of class selection to the sampleSelection package. You can get (write) access to the source code of this package on R-Forge [1]. Please let me (and Ott) know if you need any assistance. [1] http://r-forge.r-project.org/projects/sampleselection/ /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
[R] [R-pkgs] micEcon split into miscTools, micEconAids, and micEcon
The micEcon package has been split into three packages: miscTools, micEconAids, and micEcon. a) miscTools (version 0.6-0) includes miscellaneous small tools and utilities that are not related to (micro)economics, e.g. colMedians(), rowMedians(), insertCol(), insertRow(), vecli(), symMatrix(), triang(), semidefiniteness(), compPlot(), and rSquared(). The miscTools package should depend on (almost) no other packages so that it can be used in other packages without increasing the total (recursive!) dependencies of these packages too much. b) micEconAids (version 0.6-0) includes all functions and methods for demand analysis with the Almost Ideal Demand System (AIDS) of Deaton and Muellbauer (1980). Since these functions and methods comprised a large and rather specific group within the old micEcon package, it makes sense (IMHO) to separate them into a separate package. c) micEcon (version 0.6-0) includes the remaining parts of the old micEcon package, e.g. for economic analysis with Cobb-Douglas functions, translog functions, quadratic functions, CES (constant elasticity of substitution) functions, SNQ (symmetric normalized quadratic = generalized McFadden) profit functions, or linearly homogeneous non-parametric functions. I have uploaded all three packages to CRAN. They should be available for download on CRAN and its mirrors (at latest) in a few days. Feedback is highly welcome! /Arne -- Arne Henningsen http://www.arne-henningsen.name ___ R-packages mailing list r-packa...@r-project.org https://stat.ethz.ch/mailman/listinfo/r-packages __ R-help@r-project.org mailing list 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.
Re: [R] GMM estimation
Hi Zi! On 11/28/09, Yi Du abraham...@gmail.com wrote: I'd like to resort generalized methods of moment to estimate a regression model but I can't understand the help file of gmm package. If the regression model is y, the instrumental variable is x, how can I write the code? By the way, I use systemfit's 3sls with method3sls='GMM', but I can't get the suitable results. The argument method3sls of systemfit is only relevant for three-stage least-squares estimations (i.e. argument method is 3SLS). Currently, systemfit cannot estimate models by GMM. Best, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] Systemfit package
Hi Axel, On Tue, Oct 20, 2009 at 7:09 PM, Axel Leroix axel.ler...@yahoo.fr wrote: I estimated a system of equation bu using SUR method. The function summary(xx) gives me summary of estimated equation system. However, this function does not give my the value of the durbin watson statistic for each one of my equations (to chek for serial correlation). Thus, my question is is there any function in the systemfit package which permit to return the value of durbin watson statistic or should I write my own program ? The systemfit package does not provide any tools for the durbin watson test statistic (yet). So, you are invited to write this tool, whereas you could either add this via the SVN repo on R-Forge or send it to me by email :-) /Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] R formula
Hi Anwesha On Tue, Aug 18, 2009 at 10:48 AM, Anwesha Chakrabartic.anwe...@gmail.com wrote: I was trying to estimate simultaneous equation system in R using systemfit. I used the following commands library(systemfit) data(Kmenta) attach(Kmenta) eqDemand-consump~price+income eqSupply-consump~price+farmprice+trend Please replace farmprice by farmPrice fitsur-systemfit(SUR,list(demand=eqDemand, supply=eqSupply)) and got the following error messege Error in systemfit(SUR, list(demand = eqDemand, supply = eqSupply)) : argument 'formula' must be an object of class 'formula' or a list of objects of class 'formula' The error message is correct. If you do not name the arguments, the first argument must be the formula or the list of formulas. You could estimate the system either by R fitsur - systemfit( list( demand=eqDemand, supply=eqSupply ), SUR ) or by R fitsur - systemfit( method=SUR, formula=list( demand=eqDemand, supply=eqSupply ) ) Please read the documentation/manual of systemfit and the paper about systemfit published in the JSS: http://www.jstatsoft.org/v23/i04/ Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
Re: [R] SystemFit
Hi Ferdogan, Sorry for the late response. On Thu, Jul 23, 2009 at 8:29 AM, Ferdoganferdog...@hotmail.com wrote: I have two products which are substitudes. I try to fix a system as below to mydata. Demand1 = A1 -B1*Price1 + C1*Price2 Demand2 = A2 +B2*Price1 - C2*Price2 I would expect C1 B2 to be symmetric, If they are truly substitude. How can I enforce this symmetry when creating a system of equations via SystemFit ? Please take a look at http://www.jstatsoft.org/v23/i04/ I suggest that you estimate the system without the restriction first and then try to impose the restriction (and test it!). If you still do not manage to estimate the (restricted) system, please don't hesitate to send another email with a more precise/specific question to this mailing list. Best wishes, Arne -- Arne Henningsen http://www.arne-henningsen.name __ R-help@r-project.org mailing list 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.
