Home › Forums › General

Linear regression for non-overlapping segments of a wave

‹‹ Previous Topic Next Topic ››
Eduard Kas June 9, 2012 - 11:24
Posts: 40
 Joined: 2009-08-11
 Location: Netherlands

To analyse the persistence of share prices I have written a set of functions for the first time in my life. The function appears to work well and as always I am very impressed about the power of Igor Pro. Nevertheless, I need to re-write the last function (Fluctuationmagnitude), so that ultimately the regression analysis can be performed for segments of the wave instead of the entire wave, but I fail to see how this might be achieved.

What I would like to do is calculate the regression line first of all for a segment with the first 32 points, then for the subsequent 32 points, next for the third segment of 32 points, and so on, up to the 30th segment.

Per segment I would also like to calculate the Fluctuation magnitude (see variable Fluctuationmagn) and write the outcomes (a single real number for each of the 30 segments) to a wave, so that another curve fit can be performed. This time with the Fluctuation magnitude as the dependent factor and the segment number as the independent factor.

Any help is highly welcome. Many thanks in advance.

Eduard Kas 

<code>#pragma rtGlobals=1		// Use modern global access method.
 
Function PriceTrend(Startlevel)
	Variable Startlevel  	// Price of share at time = 0
 
	Variable Simulation = CreatePrices(Startlevel) 	// Fuction to simulate random walk of share 
													// prices	
	Variable Returncalc =  CalcReturn(input)		// Function to calculate  natural logs of return
												// per period	
	Variable Demeanedcalc = Demeaned(periodreturn)		// Function to calculate demeaned returns	
	Variable FFactor = Fluctuationmagnitude(squareddiff) 	// Function to calculate local trend
															// and fluctuation magnitude 
End
 
 
Function CreatePrices(Startlevel)	// Function to simulate random walk of share prices 
	Variable Startlevel	//Share price at time = 0 in euros	
 
	Make/O/N=961  input	// Create wave for share prices 
 
	Input[0] = Startlevel	// Share price at time = 0	
	Input[1,960] =Input(p-1) + gnoise(0.15)	// Share prices in subsequent period, where price is 
											// price in preceding period plus period return drawn
											// from a Gaussian (noise) distribution 
End
 
 
Function CalcReturn(input)	// Function to calculate natural logs of return per period 
	Wave Input	// Read share prices  generated by preceding function 
	String OutputName = NameOfWave(input)+"_logarithm"		// Store the name of wave plus
																// "logartihm" into the new wave 
																// OuputName 
 
	Duplicate/O input $OutputName	// Duplicate the wave with generated share prices as the 
									// wave named input_logarithm 
 
	Wave output = $OutputName		// Reference input_logarithm so that you can use it 	
	output = Ln(input)	// Calculate logarithmic values of wave 
 
	Make/O/N=960 Periodreturn
 
	Periodreturn[0]= 0 	// Assume return in first period to be zero
	Periodreturn[1,959] = output(p) - output(p-1)	// In next period logartihmic return equals logarithmic
											//  share price of this period minus logarithmic 
											// share price in previous period 
End
 
 
 
 
 
Function Demeaned(periodreturn)	// Function to calculate demeaned returns
	Wave periodreturn	// Read return values generated by previous function
 
	WaveStats/Q periodreturn	// Determine statistis of return wave
	Variable avg = V_avg	// Determine average return
 
	Make/O/N=960 Cumulative	// Create wave that will include each period's return above or
								//  below the average for all periods
	Cumulative[0] = Periodreturn[0]-V_avg	// Determine return below or above the average for 
										// the first period
	Cumulative[1,959] = Cumulative(p-1) + Periodreturn(p) - V_avg	// Determine return below or
																// above the average for all 
																// subsequent periods
End
 
 
Function Fluctuationmagnitude(cumulative)	// Function to calculate local trend and fluctuation magnitude  
	Wave cumulative	// Read cumulative demeaned returns from previous function 
 
	Make/O/N=960 Number 	// Create a new wave to show the period numbers
 
	Variable i = 0	// Initialise i
 
	do 	
 
	Number[i] = i +1 	// Add the number 1 to the level in the preceding period, so that the first
						// period carries the number 1  
 
	i += 1
 
	while(i+1<959)
 
	CurveFit line cumulative[0,959] /x=NUMBER /D // Generate a linear regression function in which 
													// the y values (demeaned cumulative 
													// returns are regressed against the x values) 
													// (the period numbers)
 
 
	Wave Fitted = fit_Cumulative	// Read fitted regression line into Wave named 'Fitted'
 
	Wave Squareddiff	// Create wave squared differences 
 
	Squareddiff = (Cumulative - Fitted)^2	// Calculate squared differences between
														//  cumulative demeaned returns and fitted 
														// regression line 
 
	WaveStats Squareddiff	// Dtermine statistics of wave squared differences 
 
	Variable Sumofsquareddiff 	// Create parameter sum of squared differences 
 
	Sumofsquareddiff = V_Sum //Determine sum
 
	Variable Fluctuationmagn 	// Create wave Fluctuation magnitude 
 
	Fluctuationmagn = Sqrt(V_Sum/960)	// Calculate fluctuation magnitude 
End<igor><code>


[ last edited June 9, 2012 - 11:42 ]
hrodstein
 June 11, 2012 - 13:23
Posts: 907
 Joined: 2007-06-21
 Location: United States

Quote:
What I would like to do is calculate the regression line first of all for a segment with the first 32 points, then for the subsequent 32 points, next for the third segment of 32 points, and so on, up to the 30th segment.

You can limit a curve fit to a subset of a wave using square brackets like this:

Make/O testWave = p + gnoise(5)
Display testWave
CurveFit line  testWave[40,80] /D 
ModifyGraph rgb(fit_testWave)=(0,0,65535), lsize(testWave)=3

I recommend changing your Fluctuationmagnitude to take startPoint, endPoint parameters passed in from a new, higher-level function. The new function would create the wave now created in Fluctuationmagnitude. It would then go into a loop calling Fluctuationmagnitude 30 times, once for each segment, passing startPoint and endPoint as parameters.

In the higher level function you can create a 30-point wave to store the outcome from each invocation of Fluctuationmagnitude. Fluctuationmagnitude would return this outcome result.

BTW, this:

	Make/O/N=960 Number 	// Create a new wave to show the period numbers
	Variable i = 0	// Initialise i
	do 	
		Number[i] = i +1 	// Add the number 1 to the level in the preceding period, so that the first
						// period carries the number 1  
	i += 1
	while(i+1<959)

can be written like this:

Make/O/N=960 Number = p+1

For details, execute:

DisplayHelpTopic "Waveform Arithmetic and Assignment"

Also, this:

Wave Fitted = fit_Cumulative	// Read fitted regression line into Wave named 'Fitted'

does not read anything. It creates a reference to fit_Cumulative.

Likewise, this does not create any wave. Rather it creates a reference to an existing wave:

Wave Squareddiff	// Create wave squared differences

A wave reference is a symbol local to a function the references a wave. For details:

DisplayHelpTopic "Wave References"
johnweeks
 June 11, 2012 - 13:26
Posts: 836
 Joined: 2007-06-29
 Location: United States

I'm not sure where you want the result to go, but here's a skeleton of what you need to do, assuming that what you want is to save the slope and intercept:

Function SegmentedLineFit(Ywave, Xwave, seglength)
	Wave Ywave, Xwave
	Variable seglength
 
	Variable nsegments = floor(numpnts(Ywave)/seglength)
	Variable i
 
	Make/O/D/N=(nsegments, 2) fitAandB		// you probably want to make this more general by allowing a customized name
	for (i = 0; i < nsegments; i += 1)
		CurveFit/Q/W=2 line, Ywave[i*seglength, (i+1)*seglength - 1] /X=Xwave
		Wave W_coef
		fitAandB[i][] = W_coef[q]
	endfor
end

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
Eduard Kas June 17, 2012 - 20:43
Posts: 40
 Joined: 2009-08-11
 Location: Netherlands

The advice of both of you has been really helpful. On that basis I created the functions below. 

#pragma rtGlobals=1		// Use modern global access method.
 
Function HurstCalculation(cumulative, Xwave, length, segmlength,nsegments)	// Create Hurst Calculation function 	
	Wave cumulative, Xwave		// Create waves 
	Variable segmlength, length, nsegments	// Create variables 
	Variable FittedLine = SegmentedLineFit(cumulative,Xwave, length, segmlength, nsegments) 	
End
 
 
Function SegmentedLineFit (cumulative,Xwave, length, segmlength, nsegments) 	// Function to calculate segmented Curve Fit 
	Wave cumulative, Xwave		// Create parameters cumulative (returns) and Xwave (= time scale)
	Variable segmlength, length, nsegments	// Create paramaters segmlength (= length of segment), length (= length of time 
											// series for price variable) and nsegments (= number of segments)
 
	Variable i 	// Create local variable i 
 
	Make/O/N=(length-1) LogWave = Ln(XWave)	// Transform time scale into natural logarithm 
 
	Variable Summedsquaredresids	// Create local variable sum of squared residuals 
 
	Make/O/N=(length-1) Squaredresids	// Create wave to which residuals for each of the successive curve fits are written down to 
	Make/O/N=(1, 1) Fluctuation	// Create wave to which square root of average residuals is written down to 
		for (i=0; i < nsegments; i += 1)	// Repeat cuve fit calculation for each equal size segment 
			CurveFit/W=2/Q line Cumulative[i*segmlength, (i+1)*segmlength-1] /X=LogWave /R	// Carry out curve 
				Wave Res_Cumulative	// Determine residuals for each curve 
				Squaredresids = Res_Cumulative^2	// Determe squared residuals for each curve fit 
		endfor
		WaveStats/Q Squaredresids 		// Determine statistics for the entire wave Squaredresids	
		Summedsquaredresids = V_avg	// Determine average value of wave
		Fluctuation = Sqrt(V_avg) 	// Take square root of the average value 
End<igor>
 
 
Eduard


[ last edited June 17, 2012 - 20:47 ]
Eduard Kas June 17, 2012 - 20:49
Posts: 40
 Joined: 2009-08-11
 Location: Netherlands

To play with the function I would like to change the Xwave depending on the length of the segment. In case of a segment length of 16, I would like to
set the first 16 points of the Xwave equal to 1 up to and including 16 and to repeat these values 60 times for an entire series length of 960 in terms of
return values (and 961 in terms of price values). (So far I created the Function SeriesCalc separately in order to test it. Later on, I would like to integrate
it with the other functions and determine the number of segments, depending on the length of a segment.)

The Function SeriesCalc does not behave like expected. If I call it by SeriesCalc(961, 16) for each of the 960 points of Xwave the value of -928 is created.
My question: what I am doing wrong.

Any help is highly welcome. Many thanks in advance.

Eduard 

#pragma rtGlobals=1		// Use modern global access method.
 
Function SeriesCalc(serieslength, segmlength)	// Function to create variable lengths of segments 
	Variable serieslength, segmlength	// Create parameters 
 
	Make/O/N=(serieslength -1) Xwave	// Create X (time) variable equal to series length (for price series)
										// minus 1
 
	Variable nsegments = floor(serieslength-1)/segmlength	// Parameter nsemgents must be an integer
 
	Variable i, j 	// Create local variables 
 
	i=0 	
	do 
		j=0	
		do 
			Xwave = j+1 - (i*segmlength)	// Set values of Xwave equal to 1 up to latest segment value 
											// that is the value of 16 if segment length is equal to 16 
			j += 1 
		while (j < segmlength)
		i += 1
	while (i < nsegments)	// Repeat Xwave values for each of the n segments 
End<igor>


[ last edited June 17, 2012 - 21:30 ]
hrodstein
 June 17, 2012 - 22:28
Posts: 907
 Joined: 2007-06-21
 Location: United States

I'm not sure that I follow what you are trying to do but I suspect that this statement is wrong:

Xwave = j+1 - (i*segmlength)

This sets every point of Xwave every time through the loop. Perhaps you mean something like:
Xwave[i*segmlength+j] = j+1

This is equivalent to a single statement:

Make /O /N=(960) xWave = mod(p,16) +1

For details on this kind of statement:

DisplayHelpTopic "Waveform Arithmetic and Assignment"

You will probably have to read this section multiple times and do some experimentation for it to sink in.


[ last edited June 17, 2012 - 22:31 ]
Eduard Kas June 19, 2012 - 23:24
Posts: 40
 Joined: 2009-08-11
 Location: Netherlands

Thank you for your suggestions. This delivers exactly what I wished to do. It is indeed a good idea to reread the text about wave references a variety of times to let it sink in and I have started to do just that.

Eduard
‹‹ Previous Topic Next Topic ››

Back to top