www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/WaveletChangepointDetectionExample.m
%% Wavelet Changepoint Detection % This example shows how to use wavelets to detect changes in the variance % of a process. Changes in variance are important because they often % indicate that something fundamental has changed about the data-generating % mechanism. % % The first example applies wavelet changepoint detection to a very old % time series -- the Nile river minima data for the years 622 to 1281 AD. % The river-level minima were measured at the Roda gauge near Cairo. % Measurements are in meters. % % Load and plot the data. load nileriverminima; years = 622:1284; figure; plot(years,nileriverminima); title('Nile River Minimum levels'); AX = gca; AX.XTick = 622:100:1222; grid on; xlabel('Year'); ylabel('Meters'); %% % Construction began on a new measuring device around 715 AD. Examining the % data prior to and after approximately 722 AD, there appears to be a % change in the variability of the data. You can use wavelets to explore % the hypothesis that the variability of the measurements has been affected % by the introduction of a new measuring device. %% % Obtain a multiresolution analysis (MRA) of the data using the Haar % wavelet. wt = modwt(nileriverminima,'haar',4); mra = modwtmra(wt,'haar'); %% % Plot the MRA and focus on the level-one and level-two details. figure; subplot(2,1,1) plot(years,mra(1,:)); title('Level 1 Details'); subplot(2,1,2) plot(years,mra(2,:)); title('Level 2 Details'); AX = gca; AX.XTick = 622:100:1222; xlabel('Years'); %% % Apply an overall change of variance test to the wavelet coefficients. for JJ = 1:5 pts_Opt = wvarchg(wt(JJ,:),2); changepoints{JJ} = pts_Opt; end cellfun(@(x) ~isempty(x),changepoints,'uni',0) %% % Determine the year corresponding to the detected change of variance. years(cell2mat(changepoints)) %% % Split the data into two segments. The first segment includes the years % 622 to 721 when the fine-scale wavelet coefficients indicate a change in % variance. The second segment contains the years 722 to 1284. % Obtain unbiased estimates of the wavelet variance by scale. tspre = nileriverminima(1:100); tspost = nileriverminima(101:end); wpre = modwt(tspre,'haar',4); wpost = modwt(tspost,'haar',4); wvarpre = modwtvar(wpre,'haar',0.95,'table') wvarpost = modwtvar(wpost,'haar',0.95,'table') %% % Compare the results. Vpre = table2array(wvarpre); Vpost = table2array(wvarpost); Vpre = Vpre(1:end-1,2:end); Vpost = Vpost(1:end-1,2:end); Vpre(:,1) = Vpre(:,2)-Vpre(:,1); Vpre(:,3) = Vpre(:,3)-Vpre(:,2); Vpost(:,1) = Vpost(:,2)-Vpost(:,1); Vpost(:,3) = Vpost(:,3)-Vpost(:,2); figure; errorbar(1:4,Vpre(:,2),Vpre(:,1),Vpre(:,3),'ko','markerfacecolor',[0 0 0]); hold on; errorbar(1.5:4.5,Vpost(:,2),Vpost(:,1),Vpost(:,3),'b^','markerfacecolor',[0 0 1]); set(gca,'xtick',1.25:4.25); set(gca,'xticklabel',{'2 Year','4 Years', '8 Years', '16 Years','32 Years'}); grid on; ylabel('Variance'); title('Wavelet Variance 622-721 and 722-1284 by Scale','fontsize',14); legend('Years 622-721','Years 722-1284','Location','NorthEast'); %% % The wavelet variance indicates a significant change in variance between % the 622-721 and 722-1284 data over scales of 2 and 4 years. %% % The above example used the Haar wavelet filter with only two coefficients % because of concern over boundary effects with the relatively small time % series (100 samples from 622-721). If your data are approximately first % or second-order difference stationary, you can substitute the biased % estimate using the 'reflection' boundary. This permits you to use a % longer wavelet filter without worrying about boundary coefficients. % Repeat the analysis using the default 'sym4' wavelet. wpre = modwt(tspre,4,'reflection'); wpost = modwt(tspost,4,'reflection'); wvarpre = modwtvar(wpre,[],[],'EstimatorType','biased','Boundary','reflection','table'); wvarpost = modwtvar(wpost,[],[],'EstimatorType','biased','Boundary','reflection','table'); %% % Plot the results. Vpre = table2array(wvarpre); Vpost = table2array(wvarpost); Vpre = Vpre(1:end-1,2:end); Vpost = Vpost(1:end-1,2:end); Vpre(:,1) = Vpre(:,2)-Vpre(:,1); Vpre(:,3) = Vpre(:,3)-Vpre(:,2); Vpost(:,1) = Vpost(:,2)-Vpost(:,1); Vpost(:,3) = Vpost(:,3)-Vpost(:,2); figure; errorbar(1:4,Vpre(:,2),Vpre(:,1),Vpre(:,3),'ko','markerfacecolor',[0 0 0]); hold on; errorbar(1.5:4.5,Vpost(:,2),Vpost(:,1),Vpost(:,3),'b^','markerfacecolor',[0 0 1]); set(gca,'xtick',1.25:4.25); set(gca,'xticklabel',{'2 Years','4 Years', '8 Years', '16 Years','32 Years'}); grid on; ylabel('Variance'); title({'Wavelet Variance 622-721 and 722-1284 by Scale'; ... 'Biased Estimate -- Reflection Boundary'},'fontsize',14); legend('622-721','722-1284','Location','NorthEast'); hold off; %% % The conclusion is reinforced. There is a significant difference in the % variance of the data over scales of 2 and 4 years, but not at longer % scales. You can conclude that something has changed about the process % variance. %% % In financial time series, you can use wavelets to detect changes in % volatility. To illustrate this, consider the quarterly chain-weighted % U.S. real GDP data for 1974Q1 to 2012Q4. The data were transformed by % first taking the natural logarithm and then calculating the % year-over-year difference. Obtain the wavelet transform (MODWT) of the % real GDP data down to level six with the 'db2' wavelet. Examine the % variance of the data and compare that to the variances by scale obtained % with the MODWT. load GDPcomponents realgdpwt = modwt(realgdp,'db2',6,'reflection'); gdpmra = modwtmra(realgdpwt,'db2','reflection'); vardata = var(realgdp,1); varwt = var(realgdpwt(:,1:numel(realgdp)),1,2); %% % In |vardata| you have the variance for the aggregate GDP time series. In % |varwt| you have the variance by scale for the MODWT. There are seven % elements in |varwt| because you obtained the MODWT down to level six % resulting in six wavelet coefficient variances and one scaling % coefficient variance. Sum the variances by scale to see that the variance % is preserved. Plot the wavelet variances by scale ignoring the scaling % coefficient variance. totalMODWTvar = sum(varwt); bar(varwt(1:end-1,:)) AX = gca; AX.XTickLabels = {'[2 4)','[4 8)','[8 16)','[16 32)','[32 64)','[64 128)'}; xlabel('Quarters') ylabel('Variance') title('Wavelet Variance by Scale') %% % Because this data is quarterly, the first scale captures variations % between two and four quarters, the second scale between four and eight, % the third between 8 and 16, and so on. % % From the MODWT and a simple bar plot, you see that cycles in the data % between 8 and 32 quarters account for the largest variance in the GDP % data. If you consider the wavelet variances at these scales, they account % for 57% of the variability in the GDP data. This means that oscillations % in the GDP over a period of 2 to 8 years account for most of the % variability seen in the time series. %% % Plot the level-one details, D1. These details capture oscillations in the % data between two and four quarters in duration. helperFinancialDataExample1(gdpmra(1,:),years,... 'Year over Year Real U.S. GDP - D1') %% % Examining the level-one details, it appears there is a reduction of % variance beginning in the 1980s. % % Test the level-one wavelet coefficients for signficant variance % changepoints. pts_Opt = wvarchg(realgdpwt(1,1:numel(realgdp)),2); years(pts_Opt) %% % There is a variance changepoint identified in 1982. This example does not % correct for the delay introduced by the 'db2' wavelet at level one. % However, that delay is only two samples so it does not appreciably affect % the results. % % To assess changes in the volatility of the GDP data pre and post 1982, % split the original data into pre- and post-changepoint series. Obtain the % wavelet transforms of the pre and post datasets. In this case, the series % are relatively short so use the Haar wavelet to minimize the number of % boundary coefficients. Compute unbiased estimates of the wavelet variance % by scale and plot the result. tspre = realgdp(1:pts_Opt); tspost = realgdp(pts_Opt+1:end); wtpre = modwt(tspre,'haar',5); wtpost = modwt(tspost,'haar',5); prevar = modwtvar(wtpre,'haar','table'); postvar = modwtvar(wtpost,'haar','table'); xlab = {'[2Q,4Q)','[4Q,8Q)','[8Q,16Q)','[16Q,32Q)','[32Q,64Q)'}; helperFinancialDataExampleVariancePlot(prevar,postvar,'table',xlab) title('Wavelet Variance By Scale'); legend('Pre 1982 Q2','Post 1982 Q2','Location','NorthWest'); %% % From the preceding plot, it appears there are significant differences % between the pre-1982Q2 and post-1982Q2 variances at scales between 2 and % 16 quarters. % % Because the time series are so short in this example, it can be useful to % use biased estimates of the variance. Biased estimates do not remove % boundary coefficients. Use a 'db2' wavelet filter with four coefficients. wtpre = modwt(tspre,'db2',5,'reflection'); wtpost = modwt(tspost,'db2',5,'reflection'); prevar = modwtvar(wtpre,'db2',0.95,'EstimatorType','biased','table'); postvar = modwtvar(wtpost,'db2',0.95,'EstimatorType','biased','table'); xlab = {'[2Q,4Q)','[4Q,8Q)','[8Q,16Q)','[16Q,32Q)','[32Q,64Q)'}; figure; helperFinancialDataExampleVariancePlot(prevar,postvar,'table',xlab) title('Wavelet Variance By Scale'); legend('Pre 1982 Q2','Post 1982 Q2','Location','NorthWest'); %% % The results confirm our original finding that there is a reduction in % volatility over scales from 2 to 16 quarters. %% % Using the wavelet transform allows you to focus on scales where the % change in volatility is localized. To see this, examine a plot of the % raw data along with the level-one wavelet details. subplot(2,1,1) helperFinancialDataExample1(realgdp,years,... 'Year over Year Real U.S. GDP -- Raw Data') subplot(2,1,2) helperFinancialDataExample1(gdpmra(1,:),years,... 'Year over Year Real U.S. GDP -- Wavelet Level 1 Details') %% % The shaded region is referred to as the "Great Moderation" signifying a % period of decreased macroeconomic volatility in the U.S. beginning in the % mid 1980s. % % Examining the aggregate data, it is not clear that there is in fact % reduced volatility in this period. However, the wavelet level-one details % uncover the change in volatility.