Re: [R] STAR models estimation with R
It is likely that here STAR means the Smooth Transition Autoregressive model which generalizes the Threshold Autoregressive (TAR) model to allow for gradual switching between regimes. The model can be estimated in R by applying maximum likelihood (ML) estimation or Monte Carlo Markov Chain (MCMC) estimation. Use Google to find related papers. Regards, Hannu Kahra -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Phineas Campbell Sent: 28. kesäkuuta 2005 15:11 To: r-help@stat.math.ethz.ch Subject: Re: [R] STAR models estimation with R The only STAR model I have come across is Smooth Threshold Autoregressive time series model, see Tong, Non Linear Time Series. I am not aware of any package that has implemented threshold models. Because statistical techniques are used widely across many different disciplines it is inevitable that naming conventions will diverge, so the same technique may have different names in different areas of study. Perhaps it should be added to the posting guide that most general name of a technique should be used, and in the case of more obscure techniques a brief reference to a standard text or paper. Hamilton, see Time Series Analysis, uses the EM algorithm to estimate such models so it should be possible to do in R. HTH Phineas Campbell -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of noomen ajmi Sent: Tuesday, June 28, 2005 8:44 AM To: r-help@stat.math.ethz.ch Subject: [R] STAR models estimation with R Hi, Can you tell me if there are an R package or code for STAR model estimation and test misspecification. If no, how i could do this. Thanks in advance Best regards AJMI Noomen Phd student TUNISIA - [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] A somewhat off the line question to a log normal distribution
Sigfried, I am not a statistician, but I have learned that according to the Central Limit Theorem (CLT) sums of random variables, regardless of their form, will tend to be normally distributed. CLT does not require the variables in the sum to come from the same underlying distribution. Ciao, Hannu Kahra Progetti Speciali Monte Paschi Asset Management SGR S.p.A. Via San Vittore, 37 IT-20123 Milano, Italia Tel.: +39 02 43828 754 Mobile: +39 333 876 1558 Fax: +39 02 43828 247 E-mail: [EMAIL PROTECTED] Web: www.mpsam.it -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of Siegfried Gonzi Sent: Thursday, December 02, 2004 10:18 AM To: [EMAIL PROTECTED] Subject: [R] A somewhat off the line question to a log normal distribution Hello: Oh yes I know it isn't so much related to R, but I gather there are a lot of statisticians reading the mailing list. My boss repeatedly tried to explain me the following. == Lets assume you have got daily measurements of a variable in natural sciences. It turned out that the aformentioned daily measurements follow a log-normal distribution when considered over the course of a year. Okay. He also tried to explain me that the monthly means (based on the daily measurements) must follow a log-normal distribution too then over the course of a year. == I somehow get his explanation. But I have measurements which are log-normal distributed when evaluated on a daily basis over the course of a year but they are close to a Gaussian distribution when considered under the light of monthly means over the course of a year. Is such a latter case feasible. And if not why. Regards, Siegfried Gonzi > > > __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] gibbs sampling for mixture of normals
Matteo, have a look at http://www.mrc-bsu.cam.ac.uk/bugs/faqs/contents.shtml and BUGS http://www.mrc-bsu.cam.ac.uk/bugs/welcome.shtml. Hannu Kahra Progetti Speciali Monte Paschi Asset Management SGR S.p.A. Via San Vittore, 37 IT-20123 Milano, Italia Tel.: +39 02 43828 754 Mobile: +39 333 876 1558 Fax: +39 02 43828 247 E-mail: [EMAIL PROTECTED] Web: www.mpsam.it -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of matteo ruggiero Sent: Thursday, November 18, 2004 7:42 PM To: [EMAIL PROTECTED] Subject: [R] gibbs sampling for mixture of normals hi i'm looking for a gibbs sampling algorithm for R for the case of mixture of K normals, and in particular for the case of bivariate normals. i'd be grateful if anyone could send its own R-routine, at least for the univariate case. thank you in advance matteo __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] R glm
In Venables & Ripley: Modern Applied Statistics with S (MASS), (4th edition), on page 184 there is a table "Families and link functions" that gives you the available links with different families. The default and the only link with the gaussian family is identity. ciao, Hannu Kahra -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of Shuangge Ma Sent: Thursday, September 23, 2004 6:26 PM To: [EMAIL PROTECTED] Subject: [R] R glm Hello: would you please help me with the following glm question? for the R function glm, what I understand is: once you specify the "family", then the link function is fixed. My question is: is it possible I use, for example, "log" link function, but the estimation approach for the guassian family? Thanks, Shuangge Ma, Ph.D. * CHSCC, Department of Biostatistics * * University of Washington * * Building 29, Suite 310, 6200 NE 74th ST. * * Seattle, WA 98115* * Tel: 206-685-7123 Fax: 206-616-4075 * __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] block statistics with POSIX classes
I have followed Gabor's instructions: > aggregate(list(y=y), list(dp$year), mean)$y # returns NULL since y > is a time series NULL > aggregate(list(y=as.vector(y)), list(dp$year), mean)$y# returns annual means [1] 0.0077656696 0.0224050294 0.001898 0.0240550925 -0.0084085867 [6] -0.0170950194 -0.0355641251 0.0065873997 0.0008253111 > aggregate(list(y=y), list(dp$year), mean) # returns the same as > the previous one Group.1 Series.1 1 96 0.0077656696 2 97 0.0224050294 3 98 0.001898 4 99 0.0240550925 5 100 -0.0084085867 6 101 -0.0170950194 7 102 -0.0355641251 8 103 0.0065873997 9 104 0.0008253111 Gabor's second suggestion returns different results: > aggregate(ts(y, start=c(dp$year[1],dp$mon[1]+1), freq = 12), nfreq=1, mean) Time Series: Start = 96.3 End = 103. Frequency = 1 Series 1 [1,] 0.016120895 [2,] 0.024257131 [3,] 0.007526997 [4,] 0.017466118 [5,] -0.016024846 [6,] -0.017145159 [7,] -0.036047765 [8,] 0.014198501 > aggregate(y, 1, mean) # verifies the result above Time Series: Start = 1996.333 End = 2003.333 Frequency = 1 Series 1 [1,] 0.016120895 [2,] 0.024257131 [3,] 0.007526997 [4,] 0.017466118 [5,] -0.016024846 [6,] -0.017145159 [7,] -0.036047765 [8,] 0.014198501 The data is from 1996:5 to 2004:8. The difference of the results must depend on the fact that the beginning of the data is not January and the end is not December? The first two solutions give nine annual means while the last two give only eight means. The block size in the last two must be 12 months, as is said in ?aggregate, instead of a calender year that I am looking for. Gabor's first suggestion solved my problem. Thank you, Hannu -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of Gabor Grothendieck Sent: Thursday, September 23, 2004 3:52 PM To: [EMAIL PROTECTED] Subject: Re: [R] block statistics with POSIX classes I am not sure that I understand what you are looking for since you did not provide any sample data with corresponding output to clarify your question. Here is another guess. If y is just a numeric vector with monthly data and dp is a POSIXlt vector of the same length then: aggregate(list(y=y), list(dp$year), mean)$y will perform aggregation, as will aggregate(ts(y, start=c(d$year[1],d$mon[1]+1), freq = 12), nfreq=1, mean) which converts y to ts and then performs the aggregation. The first one will work even if y is irregular while the second one assumes that y must be monthly. The second one returns a ts object. By the way, I had a look at gev's source and it seems that despite the documentation it does not use POSIXct anywhere internally. If the block is "year" or other character value then it simply assumes that whatever datetime class is used has an as.POSIXlt method. If your dates were POSIXct rather than POSIXlt then it would be important to ensure that whatever timezone is assumed (which I did not check) in the conversion is the one you are using. You could use character dates or Date class to avoid this problem. Since you appear to be using POSIXlt datetimes from the beginning I think you should be ok. Kahra Hannu mpsgr.it> writes: : : Thank you Petr and Gabor for the answers. : : They did not, however, solve my original problem. When I have a monthly time series y with a POSIX date : variable dp, the most obvious way to compute e.g. the annual means is to use aggregate(y, 1, mean) that : works with time series. However, I got stuck with the idea of using the 'by' argument as by = dp$year. But in : that case y has to be a data.frame. The easiest way must be the best way. : : Regards, : Hannu : : -Original Message- : From: r-help-bounces stat.math.ethz.ch : [mailto:r-help-bounces stat.math.ethz.ch]On Behalf Of Gabor Grothendieck : Sent: Thursday, September 23, 2004 12:56 PM : To: r-help stat.math.ethz.ch : Subject: Re: [R] block statistics with POSIX classes : : : Kahra Hannu mpsgr.it> writes: : : : : : I have a monthly price index series x, the related return series y = diff (log : (x)) and a POSIXlt date-time : : variable dp. I would like to apply annual blocks to compute for example : annual block maxima and mean of y. : : : : When studying the POSIX classes, in the first stage of the learning curve, I : computed the maximum drawdown : : of x: : : > mdd <- maxdrawdown(x) : : > max.dd <- mdd$maxdrawdown : : > from <- as.character(dp[mdd$from]) : : > to <- as.character(dp[mdd$to]) : : > from; to : : [1] "2000-08-31" : : [1] "2003-03-31" : : that gives me the POSIX dates of the start and end of the period and : suggests that I have done something correctly. : : : : Two questions: : : (1) how to imple
RE: [R] block statistics with POSIX classes
Thank you Petr and Gabor for the answers. They did not, however, solve my original problem. When I have a monthly time series y with a POSIX date variable dp, the most obvious way to compute e.g. the annual means is to use aggregate(y, 1, mean) that works with time series. However, I got stuck with the idea of using the 'by' argument as by = dp$year. But in that case y has to be a data.frame. The easiest way must be the best way. Regards, Hannu -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of Gabor Grothendieck Sent: Thursday, September 23, 2004 12:56 PM To: [EMAIL PROTECTED] Subject: Re: [R] block statistics with POSIX classes Kahra Hannu mpsgr.it> writes: : : I have a monthly price index series x, the related return series y = diff(log (x)) and a POSIXlt date-time : variable dp. I would like to apply annual blocks to compute for example annual block maxima and mean of y. : : When studying the POSIX classes, in the first stage of the learning curve, I computed the maximum drawdown : of x: : > mdd <- maxdrawdown(x) : > max.dd <- mdd$maxdrawdown : > from <- as.character(dp[mdd$from]) : > to <- as.character(dp[mdd$to]) : > from; to : [1] "2000-08-31" : [1] "2003-03-31" : that gives me the POSIX dates of the start and end of the period and suggests that I have done something correctly. : : Two questions: : (1) how to implement annual blocks and compute e.g. annual max, min and mean of y (each year's max, min, mean)? : (2) how to apply POSIX variables with the 'block' argument in gev in the evir package? : : The S+FinMetrics function aggregateSeries does the job in that module; but I do not know, how handle it in R. : My guess is that (1) is done by using the function aggregate, but how to define the 'by' argument with POSIX variables? 1. To create a ts monthly time series you specify the first month and a frequency of 12 like this. z.m <- ts(rep(1:6,4), start = c(2000,1), freq = 12) z.m # Annual aggregate is done using aggregate.ts with nfreq = 1 z.y <- aggregate(z.m, nfreq = 1, max) z.y # To create a POSIXct series of times using seq # (This will use GMT. Use tz="" arg to ISOdate if you want current tz.) z.y.times <- seq(ISOdate(2000,1,1), length = length(z.y), by = "year") z.y.times 2. Have not used evir but looking at ?gev it seems you can use block = 12 if you have monthly data and want the blocks to be successive 12 month periods or you can add a POSIXct times attribute to your data as below (also see comment re tz above) and then use block = "year" in your gev call. attr(z.m, "times") <- seq(ISOdate(2000,1,1), length=length(z.m), by="month") str(z.m) # display z.m along with attribute info __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] block statistics with POSIX classes
I have a monthly price index series x, the related return series y = diff(log(x)) and a POSIXlt date-time variable dp. I would like to apply annual blocks to compute for example annual block maxima and mean of y. When studying the POSIX classes, in the first stage of the learning curve, I computed the maximum drawdown of x: > mdd <- maxdrawdown(x) > max.dd <- mdd$maxdrawdown > from <- as.character(dp[mdd$from]) > to <- as.character(dp[mdd$to]) > from; to [1] "2000-08-31" [1] "2003-03-31" that gives me the POSIX dates of the start and end of the period and suggests that I have done something correctly. Two questions: (1) how to implement annual blocks and compute e.g. annual max, min and mean of y (each year's max, min, mean)? (2) how to apply POSIX variables with the 'block' argument in gev in the evir package? The S+FinMetrics function aggregateSeries does the job in that module; but I do not know, how handle it in R. My guess is that (1) is done by using the function aggregate, but how to define the 'by' argument with POSIX variables? Thanks! Hannu Kahra Progetti Speciali Monte Paschi Asset Management SGR S.p.A. Via San Vittore, 37 IT-20123 Milano, Italia Tel.: +39 02 43828 754 Mobile: +39 333 876 1558 Fax: +39 02 43828 247 E-mail: [EMAIL PROTECTED] Web: www.mpsam.it __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] unable to install dse package
>Error in file(file, "r") : unable to open connection >In addition: Warning message: >cannot open file 'a:/project.txt' >please help help(read.table) Hannu Kahra - [[alternative HTML version deleted]] __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Estimating parameters for a bimodal distribution
George,
Venables & Ripley: Modern Applied Statistics with S, Springer 2002, Chapter 16 (An
example: fitting a mixture model) may be helpful.
Hannu Kahra
-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] Behalf Of George W.
Gilchrist
Sent: Thursday, September 16, 2004 5:00 PM
To: [EMAIL PROTECTED]
Subject: [R] Estimating parameters for a bimodal distribution
For several years, I have been using Splus to analyze an ongoing series of
datasets that have a bimodal distribution. I have used the following
functions, in particular the ms() function, to estimate the parameters: two
means, two standard deviations, and one proportion. Here is the code I've
been using in S:
btmp.bi <- function(vec, p, m1, m2, sd1, sd2)
{
(exp(p)/(1+exp(p)))*dnorm(vec,mean=m1,sd=abs(sd1))+
(1-(exp(p)/(1+exp(p*dnorm(vec,mean=m2,sd=abs(sd2))
}
btmp11 <- ms( ~ - sum(log((btmp.bi(btmp1$Temp, p, m1, m2,
s1, s2, start = list(p = 0.4, m1 = 38, m2 = 40, s1
= 1, s2 = 1), control = list(maxiter = 200))
I have looked in the archives and tried various alternatives, especially
optim(), but so far have had nothing but frustration. I've been running this
in a semi-automated program on several hundred datasets a few times each
month for several years now. I would like to figure out how to move this to
R. Thank you for any help you might offer.
==
George W. GilchristEmail #1: [EMAIL PROTECTED]
Department of Biology, Box 8795 Email #2: [EMAIL PROTECTED]
College of William & MaryPhone: (757) 221-7751
Williamsburg, VA 23187-8795Fax: (757) 221-6483
http://gwgilc.people.wm.edu/
__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Bessel function
>Currently, I'm implementing the generalized hyperbolic distribution into >Splus. Unfortunately the Bessel function is not implemented in Splus. In >R the Bessel function does exist but it is an internal function and I'm >not able to look at the code. >Is there any possibility to see the code of the Bessel function in R or >does anybody has an implementation of the Bessel function in Splus? Have a look at Press-Teukolsky-Vetterling-Flannery: Numerical Recipes in C, Cambridge University Press, Chapter 6 (2nd edition). I guess you are looking for modified Bessel functions (see, 6.6). Hannu Kahra Progetti Speciali Monte Paschi Asset Management SGR S.p.A. Via San Vittore, 37 IT-20123 Milano, Italia Tel.: +39 02 43828 754 Mobile: +39 333 876 1558 Fax: +39 02 43828 247 E-mail: [EMAIL PROTECTED] Web: www.mpsam.it __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] linear constraint optim with bounds/reparametrization
>From Spencer Graves: >However, for an equality constraint, I've had good luck by with an objective function >that adds something like the >following to my objective function: constraintViolationPenalty*(A%*%theta-c)^2, where >"constraintViolationPenalty" is >passed via "..." in a call to optim. I applied Spencer's suggestion to a set of eight different constrained portfolio optimization problems. It seems to give a usable practice to solve the portfolio problem, when the QP optimizer is not applicable. After all, practical portfolio management is more an art than a science. >I may first run optim with a modest value for constraintViolationPenalty then restart >it with the output of the >initial run as starting values and with a larger value for >constraintViolationPenalty. I wrote a loop that starts with a small value for the penalty and stops when the change of the function value, when increasing the penalty, is less than epsilon. I found that epsilon = 1e-06 provides a reasonable accuracy with respect to computational time. Spencer, many thanks for your suggestion. Hannu Kahra __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] Help with generating data from a 'not quite' Normal distriburtion
Consider using the HyperbolicDist package. With the package you can both fit the hyperbolic distribution to your data and generate random numbers from the distribution. Hyperbolic distribution/s provide/s good fit to financial returns that commonly exhibit high peaks and heavy tails. Hannu Kahra -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] Behalf Of Crabb, David Sent: Thursday, August 12, 2004 10:44 AM To: [EMAIL PROTECTED] Subject: [R] Help with generating data from a 'not quite' Normal distriburtion I would be very grateful for any help from members of this list for what might be a simple problem... We are trying to simulate the behaviour of a clinical measurement in a series of computer experiments. This is simple enough to do in R if we assume the measurements to be Gaussian, but their empirical distribution has a much higher peak at the mean and the distribution has much longer tails. (The distribution is quite symmetrical) Can anyone suggest any distributions I could fit to this data, and better still how I can then generate random data from this 'distribution' using R? --- Dr. David Crabb School of Science, The Nottingham Trent University, Clifton Campus, Nottingham. NG11 8NS Tel: 0115 848 3275 Fax: 0115 848 6690 __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html __ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
RE: [R] linear constraint optim with bounds/reparametrization
> from Ingmar Visser:
>I would like to optimize a (log-)likelihood function subject to a number of
>linear constraints between parameters. These constraints are equality
>constraints of the form A%*%theta=c, ie (1,1) %*% 0.8,0.2)^t = 1 meaning
>that these parameters should sum to one. Moreover, there are bounds on the
>individual parameters, in most cases that I am considering parameters are
>bound between zero and one because they are probabilities.
>My problems/questions are the following:
>1) constrOptim does not work in this case because it only fits inequality
>constraints, ie A%*%theta > = c
---
I was struggling with the same problem a few weeks ago in the portfolio optimization
context. You can impose equality constraints by using inequality constraints >= and <=
simultaneously. See the example bellow.
>As a result when providing starting values constrOptim exits with an error
>saying that the initial point is not feasible, which it isn't because it is
>not in the interior of the constraint space.
In constrOptim the feasible region is defined by ui%*%theta-ci >=0. My first attempt
to obtain feasible starting values was based on solving for theta from ui%*%theta =
ci. Some of the items in the right hand side of the feasible region are, however, very
often very small negative numbers, and hence, theta is not feasible. Next, I started
to increase ci by epsilon ("a slack variable") and checked if the result was feasible.
The code is in the example bellow. If ui is not a square matrix, theta is obtained by
simulation. It is helpfull to know the upper and lower limits of the theta vector. In
my case they are often [0,1]. Usually only 2-3 simulations are required.
In the example, the weights (parameters) have limits [0,1] such that their sum is
limited to unity, as in your case. See, how Amat and b0 are defined.
V <- matrix(c(0.12,0,0,0,0.21,0,0,0,0.28),3,3) # variances
Cor <- matrix(c(1,0.25,0.25,0.25,1,0.45,0.05,0.45,1),3,3) # correlations
sigma <- t(V)%*%Cor%*%V # covariance matrix
ER <- c(0.07,0.12,0.18) # expected returns
Dmat <- sigma # adopted from solve.QP
dvec <- rep(0,3)#" " "
k <- 3 # number of assets
reslow <- rep(0,k) # lower limits for
portfolio weights
reshigh <- rep(1,k) # upper limits for
portfolio weights
pm <- 0.10 # target return
rfree <- 0.05 # risk-free return
risk.aversion <- 2.5# risk aversion parameter
### RISKLESS = FALSE; SHORTS = TRUE ; CONSTRAINTS = TRUE
a1 <- rep(1, k)
a2 <- as.vector(ER)+rfree
a3 <- matrix(0,k,k)
diag(a3) <- 1
Amat <- t(rbind(a1, a2, a3, -a3))
b0 <- c(1,pm,reslow, -reshigh)
objectf.mean <- function(x) # object
function
{
return(sum(t(x)%*%ER-1/2*risk.aversion*t(x)%*%sigma%*%x)-(t(x)%*%ER-pm))
}
# Getting starting values for constrOptim
ui <- t(Amat)
ci <- b0
if (dim(ui)[1] == dim(ui)[2])
{
test <- F
cieps <- ci
while(test==F)
{
theta <- qr.solve(ui,cieps)
foo <- (ui%*%theta-ci) # check if the initial values are in
the feasible area
test <- all(foo > 0)
if(test==T) initial <- theta# initial values for portfolio weights
cieps <- cieps+0.1
}
}
if (dim(ui)[1] != dim(ui)[2]) # if Amat is not square, simulate the
starting values
{
test <- F
i <- 1
while(test==F)
{
theta <- runif(k, min = 0, max = 1)
foo <- (ui%*%theta-ci)
test <- all(foo > 0)# and check that theta is feasible
if (test == T) initial <- theta
i <- i+1
}
}
initial
res <- constrOptim(initial, objectf.mean, NULL, ui=t(Amat), ci=b0, control =
list(fnscale=-1))
res$par # portfolio weights (parameters)
y <- t(res$par)%*%ER# return on the portfolio
y
I hope this helps.
Hannu Kahra
Progetti Speciali
Monte Paschi Asset Management SGR S.p.A.
Via San Vittore, 37 - 20123 Milano, Italia
Tel.: +39 02 43828 754
Mobile: +39 333 876 1558
Fax: +39 02 43828 247
E-mail: [EMAIL PROTECTED]
Web: www.mpsam.it
"Le informazioni trasmesse sono da intendersi per esclusivo uso del destinatario e
possono contenere informazio
