@@ -87,6 +87,7 @@ impl SamplingExt<f64> for AGARCH<f64> {
8787 x
8888 }
8989
90+ #[ cfg( feature = "simd" ) ]
9091 fn sample_simd ( & self ) -> Array1 < f64 > {
9192 use crate :: stats:: distr:: normal:: SimdNormal ;
9293
@@ -135,6 +136,60 @@ impl SamplingExt<f64> for AGARCH<f64> {
135136 x
136137 }
137138
139+ fn n ( & self ) -> usize {
140+ #[ cfg( not( feature = "simd" ) ) ]
141+ self . n
142+ }
143+
144+ #[ cfg( feature = "simd" ) ]
145+ fn sample_simd ( & self ) -> Array1 < f32 > {
146+ use crate :: stats:: distr:: normal:: SimdNormal ;
147+
148+ let p = self . alpha . len ( ) ;
149+ let q = self . beta . len ( ) ;
150+
151+ // Generate white noise
152+ let z = Array1 :: random ( self . n , SimdNormal :: new ( 0.0 , 1.0 ) ) ;
153+
154+ // Arrays for X_t and sigma_t^2
155+ let mut x = Array1 :: < f32 > :: zeros ( self . n ) ;
156+ let mut sigma2 = Array1 :: < f32 > :: zeros ( self . n ) ;
157+
158+ // Summation for unconditional variance init
159+ let sum_alpha: f32 = self . alpha . iter ( ) . sum ( ) ;
160+ let sum_delta_half: f32 = self . delta . iter ( ) . sum :: < f32 > ( ) * 0.5 ;
161+ let sum_beta: f32 = self . beta . iter ( ) . sum ( ) ;
162+ let denom = ( 1.0 - sum_alpha - sum_delta_half - sum_beta) . max ( 1e-8 ) ;
163+
164+ for t in 0 ..self . n {
165+ if t == 0 {
166+ sigma2[ t] = self . omega / denom;
167+ } else {
168+ let mut var_t = self . omega ;
169+ // p-lag terms
170+ for i in 1 ..=p {
171+ if t >= i {
172+ let x_lag = x[ t - i] ;
173+ let indicator = if x_lag < 0.0 { 1.0 } else { 0.0 } ;
174+
175+ var_t +=
176+ self . alpha [ i - 1 ] * x_lag. powi ( 2 ) + self . delta [ i - 1 ] * x_lag. powi ( 2 ) * indicator;
177+ }
178+ }
179+ // q-lag terms
180+ for j in 1 ..=q {
181+ if t >= j {
182+ var_t += self . beta [ j - 1 ] * sigma2[ t - j] ;
183+ }
184+ }
185+ sigma2[ t] = var_t;
186+ }
187+ x[ t] = sigma2[ t] . sqrt ( ) * z[ t] ;
188+ }
189+
190+ x
191+ }
192+
138193 fn n ( & self ) -> usize {
139194 self . n
140195 }
0 commit comments