1

I was trying to do a linear model analysis where the parameters are constrained (sum to 1 and non-negative).But I found that is not obvious how to apply AIC function(or others) to the parameters selection as I wanted to use optimization to solve the linear problem.

Any idea or other suggestions?

thanks,

Wei

Wei
  • 13
  • 2
  • You may want to look at the LASSO instead of the AIC. – gung - Reinstate Monica Feb 05 '16 at 09:17
  • Given you can calculate the log-likelihood as well as count the parameters of the optimized model what stops you from using AIC? – usεr11852 Feb 05 '16 at 09:17
  • @usεr11852 I's new to constrained LM, so I used solve.QP to do optimization, I did not calculate the log-likelihood(if I do, how do I do it). I was thinking of finding sth like "stepAIC" that can help me choosing the model automatically. What do you say? – Wei Feb 05 '16 at 09:43
  • @gung Not familiar with LASSO, pardon my ignorance on the topic... Would you please give me a more detailed explanation? – Wei Feb 05 '16 at 09:49
  • @usεr11852 is the $Lagrangian from the output of the optimization Log-likelihood? – Wei Feb 05 '16 at 11:22
  • @Wei: Depends what you optimize against. Maybe you optimize against the deviance for example (which is a scaled log-likelihood), maybe the RSS I cannot tell. Please add you optimization procedure and/or more information about how you fit this model so we can help more substantially. – usεr11852 Feb 05 '16 at 17:09
  • @usεr11852 It has to do with the return based analysis: Y(i)= a+W(i)*X(i), where sum(W(i))=1 and Wi >=0 , i= 1,2,...40 This is the regression I wanted to do. Help please! – Wei Feb 05 '16 at 17:22
  • solve.QP(Dmat, dvec, Amat, bvec, meq=0, factorized=FALSE) Dmat=Cov(X), dvec=Y^T*X, C – Wei Feb 05 '16 at 17:29
  • Please edit your question accordingly outlining what you tried to do and why. I think that maybe a solution is possible but I cannot tell from the current problem formulation. – usεr11852 Feb 05 '16 at 18:04
  • (BTW: I think this has the potential of being an interesting question; just at this point there too many loose ends with it.) – usεr11852 Feb 05 '16 at 18:06
  • @usεr11852 The following equation is what Im trying to [Model](https://en.wikipedia.org/wiki/Returns-based_style_analysis) $$R_t^m=\alpha+\sum\limits_{i=1}^I\beta^iR_t^i+\epsilon_t $$ where: -$R_t^m $is the time stream of historical manager returns, -$R_t^i$ is a set of time streams of market indices or factors -$I $is the number of indices or factors used in analysis, -$ \alpha $is the intercept of the regression equation, often interpreted as manager skill, -$\epsilon_t $is the error, to be minimized using ordinary least squares regression. – Wei Feb 05 '16 at 22:54
  • @usεr11852 $ \beta^i $ >=0 and $ sum(\beta^i) $ =1. In my problem, $I$ is up to 40, and the goal is to reduce the number of factors $ R^i $ to less than 15 respecting the constraints and using AIC – Wei Feb 05 '16 at 23:03
  • @usεr11852 it is indeed a very interesting application, hopefully it can be sorted out – Wei Feb 05 '16 at 23:04
  • @Wei: Please see my answer below. I think it address what you wanted which in general not very common. – usεr11852 Feb 06 '16 at 09:26

1 Answers1

1

OK, one can try the following procedure to get a AIC-like test.

  1. Solve the quadratic programming problem defined and get $x_{opt}$. Say k is the length of $x_{opt}$.
  2. Use $x_{opt}$ in a regular linear model; get fitted values and then calculate that model's residuals.
  3. Fit a Gaussian distribution to that residual vector. Usual implementation use an E-M algorithm but we will not expand on this here. Computationally one can use something like resDist = MASS::fitdistr(residVector, 'normal') in R.
  4. Use that Gaussian distribution to get the log-likelihood of the distribution object. Mathematically we would calculate the Gaussian log-likelihood as: $ \log(L(\theta)) =-\frac{|D|}{2}\log(2*\pi) -\frac{1}{2} \log(|K|) -\frac{1}{2}(y-\hat{y})^T K^{-1} (y-\hat{y})$, $K$ being the covariance structure of your model, $|D|$ being the number of points in your datasets, $\hat{y}$ the mean response and obviously $y$ being your dependent variable. Computationally one can use something like: myLL = logLik(resDist) in R.
  5. Use the number of parameters k from step 1 and the log-likelihood myLL from step 4 to get your AIC. (Or even better AICc so you somewhat correct for your finite sample.)

Caveats:

  1. Between steps 2-3 we make the leap of faith that the residuals of this model are normally distributed. Are they? Maybe they are? You need to check that. Usually deviations from normality are not catastrophic but given the constraints you have you might be very far off.
  2. The AIC you get is an upper bound for the AIC of that model. This is because your estimates for the parameters in $x_{opt}$ are not true maximum likelihood estimates. They would be only if your constrained and unconstrained solutions were equal. The fact that it is an upper bound means that it will only be good as an indicator. It is not a true AIC so when you report you should keep that in mind.
  3. On first instance do not bootstrap this model. Your $x_{opt}$ solution potentially lies close to (or even onto) your parameters' space boundary. As a consequence your bootstrap might be biased.
  4. Stepwise regression (because essentially this is what you want to do) is known to cause data-dredging which is a bad thing. It can lead to false interpretation of the results. Take any variable selection procedure with a grain of salt. This is issue is extesnsively discussed in Cross-Validated; eg. see the following two great threads for starters here and here.
usεr11852
  • 33,608
  • 2
  • 75
  • 117
  • Thank you for you explanation and solution, Ive learned a lot on this topic! Let me check if it works then I will let you know the results. Btw, you mantioned above that Stepwise was know to cause data dredging, what other criteria would you use to do model choice? – Wei Feb 08 '16 at 12:39
  • I am glad I could help; if this is post is useful for you or answers your original question please consider upvoting it or accepting it as an answer. The obvious thing to try is cross-validation first of all. The threads linked provide more alternatives. – usεr11852 Feb 08 '16 at 20:05